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Preface 


The following chapters contain core material supported by pen and paper exercises together with 
computer-based exercises where appropriate. In addition there are web links to: 


e worked solutions, 

e computer codes, 

e audio-visual presentations, 
e case studies, 

e further reading. 


Codes are written using Scilab (a Matlab clone, downloadable for free from http://www.scilab.org/) and 
also Matlab. 


The emphasis of this book is on the practical: students are encouraged to experiment with different input 
parameters and investigate outputs in the computer-based exercises. Theory is reduced to a necessary 
minimum and provided in appendices. Web links are found on the following web page: 


http://www2.docm.mmu.ac.uk/STAFF/C.Mingham/ 


This book is intended for final year undergraduates who have knowledge of Calculus and introductory 
level computer programming. 
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1. Introduction 


This book provides an introduction to the finite difference method (FDM) for solving partial differential 
equations (PDEs). In addition to specific FDM details, general concepts such as stability, boundary 
conditions, verification, validation and grid independence are presented which are important for anyone 
wishing to solve PDEs by using other numerical methods and/or commercial software packages. Material 
is presented in order of increasing complexity and supplementary theory is included in appendices. 


1.1 Partial Differential Equations 


The following equation is an example of a PDE: 
a(t, x,y) U, + b(t,x,y)U,+c(t,x,y)U,, =f(t.x,y) (1.1) 


where, 


e t,x, y are the independent variables (often time and space) 
e a,b,c and fare known functions of the independent variables, 
e Us the dependent variable and is an unknown function of the independent variables. 


oU oU OU 
¢ partial derivatives are denoted by subscripts: U,\=—,U,=—, U,,==>5 
ot Ox oy 


The order of a PDE is the order of its highest derivative. A PDE is /inear if U and all its partial derivatives 


etc. 


occur to the first power only and there are no products involving more than one of these terms. (1.1) is 
second order and linear. The dimension of a PDE is the number of independent spatial variables it 
contains. (1.1) is 2D if x and y are spatial variables. 


1.2 Solution to a Partial Differential Equation 


Solving a PDE means finding the unknown function U. An analytical (i.e. exact) solution of a PDE is a 
function that satisfies the PDE and also satisfies any boundary and/or initial conditions given with the PDE 
(more about these later). Most PDEs of interest do not have analytical solutions so a numerical procedure must 
be used to find an approximate solution. The approximation is made at discrete values of the independent 
variables and the approximation scheme is implemented via a computer program. The FDM replaces all partial 
derivatives and other terms in the PDE by approximations. After some manipulation, a finite difference scheme 
(FDS) is created from which the approximate solution is obtained. The FDM depends fundamentally on 
Taylor’s beautiful theorem (circa 1712!) which is stated in the next chapter. 
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1.3 PDE Models 


PDEs describe many of the fundamental natural laws (e.g. conservation of mass) so describe a wide range 
of physical phenomena. Examples include Laplace’s equation for steady state heat conduction, the 
advection-diffusion equation for pollutant transport, Maxwell’s equations for electromagnetic waves, the 
Navier—Stokes equations for fluid flow and many, many more. The authors’ main interest is in solving 
PDEs for fluid flow problems and details, including pictures and animations, can be found at: 
http://www.docm.mmu.ac.uk/emmfa/ 


1.4 Classification of PDEs 


Second order linear PDEs can be formally classified into 3 generic types: elliptic, parabolic and 
hyperbolic. The simplest examples are: 


a) Elliptic: e.g.U,, +U,, =f(x, y). 


This is Poisson’s equation or Laplace’s equation (when f(x,y) =0) which may be used to model the steady 
state temperature distribution in a plate or incompressible potential flow. Notice there is no time derivative. 


b) Parabolic: e.g. U, = kU,,. 


This is the 1D diffusion equation and can be used to model the time-dependent temperature distribution 
along a heated 1D bar. 


c) Hyperbolic: e.g. U, = c°U 


This is the wave equation and may be used to model a vibrating guitar string or 1D supersonic flow. 

d) U, =-cU,. 

This first order PDE is called the advection equation. Solutions of d) also satisfy c). 

e) U,+cU, =kU,,. 

This is the advection-diffusion equation and may be used to model transport of a pollutant in a river. The 
coefficients k, c in the above PDEs quantify material properties that relate to the problem being solved e.g. 


k could be the coefficient of thermal conductivity in the case of a heated bar, or 1D diffusion coefficient in 
the case of pollutant transport; c is a wave speed, usually, in fluid flow, the speed of sound. 
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1.4.1 Initial and Boundary Conditions 


PDEs require proper initial conditions (ICs) and boundary conditions (BCs) in order to define what is 
known as a well-posed problem. If too many conditions are specified then there will be no solution; if too 
few conditions are specified the solution will not be unique. If the ICs/BCs are specified in the wrong 
place or at the wrong time then the solution will not depend smoothly on the ICs/BCs and small errors in 
the ICs/BCs will bring about large changes in the solution. This is referred to as an ill-posed problem. The 
PDEs encountered in practice are often non-linear and multi-dimensional and cannot be reduced to the 
simple so-called canonical forms of a) - e). However, we need to understand the properties of the solution 
to these simple model PDEs before attempting to solve more complicated PDEs. 


A second order elliptic PDE such as a) requires a boundary condition on U at each point on the boundary. 
Thus, these are called Boundary Value (BV) problems. The BC may be a value of U on the boundary or 
the value of its derivative (see Appendix B). Linear parabolic equations such a b) require ICs at the initial 
start time (usually t=0) and one BC at each end-point of the spatial domain (e.g. at the ends of the heated 
bar). Technically linear hyperbolic equations such as d) require ICs and as many BCs as there are inward- 
pointing characteristics (this is an advanced topic which we will not cover) which depend on the sign of 
wave speed c, thus: 


If c>0, we need ICs: U(0,x) = f(x) and BCs: U(t,0) = g(t); 
If c<0, we need ICs: U(0,x) = f(x) but no BCs. 
These are called Initial Boundary Value Problems (IBV) problems. 
1.4.2 Domain of Dependence 


The differences between the types of PDEs can be illustrated by sketching their respective domains of 
dependence. So for example, in the hyperbolic case d), point P (Xo, to) in Figure 1.1 can only be influenced 
by points lying within the region bounded by the two characteristics x+ct = const and x-ct = const and 

t < to. This region is called the domain of dependence. In turn, point P can influence points at later times 
lying within its zone of influence. In the parabolic case, shown in Figure 1.2 information travels 
downstream (or forward in time) only and so the domain of dependence of point P (Xo, to) in this case is 
the region t < to and the zone of influence is all points for which t > to, 
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Zone of 
influence 


X+ct = const 


P (Xo, to) 


x-ct = const 
Domain of 
dependence 


IC | x 
Figure 1.1 Domain of dependence: hyperbolic case. 


Zone of 
influence 


BC BC 
P (Xo, to) 


Domain of 
dependence 


IC 
Figure 1.2 Domain of dependence: parabolic case. 
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In the elliptic case, corresponding to subsonic flow (Figure 1.3), information travels in all directions at 
infinite speed so the solution at point P (xo, to) influences all points within the domain and vice versa. 


BC 


P (Xo, to) 


BC Domain of Zone of 


dependence influence 


BC xX 
Figure 1.3 Domain of dependence: elliptic case. 


Notice in this case that the whole region bounded by the BCs is both a domain of dependence and zone of 


influence. 


6007 'A'G Swewsds Wy] 222] @ 


iB \ 
3 It’s only an 


opportunity if 
you act on it 


IKEA.SE/STUDENT 


on 11S fos a 


Download free ebooks at bookboon.com 


14 


Introductory Finite Difference Methods for PDEs Introduction 


The type of PDE fundamentally influences the choice of solution strategy. Time dependent hyperbolic 
problems and parabolic problems illustrated by Figures 1.1 and 1.2 are solved numerically by time- 
marching methods which involves, as its name suggests, obtaining the numerical solution at a later time 
from that at an earlier time starting from given ICs. 


Elliptic problems, as illustrated in Figure 1.3 are solved numerically by so-called relaxation methods. 


1.5 Discrete Notation 


We will use upper case U to denote the analytic (exact) solution of the PDE and lower case u to denote the 
numerical (approximate) solution. Subscripts will denote discrete points in space and superscripts discrete 


levels in time. e.g. u’ ; denotes the numerical solution at grid point (i, j) in a 2D region at time level n. 


1.6 Checking Results 


Before applying a numerical scheme to real life situations modelled by PDEs there are two important steps 
that should always be undertaken. 


1.6.1 Verification 


The computer program implementing the scheme must be verified. This is a check to see if the program is 
doing what it is supposed to do. Comparing results from pen and paper calculations at a small number of 
points to equivalent computer output is a way to (partially) verify a program. Give or take a small amount 
of rounding error the numbers should be the same. Another way to verify the program is to find an exact 
solution to the PDE for a simpler problem (if one exists) and compare numerical and exact results. 
Complete program verification involves testing that all branches, program elements and statements are 
executed and produce the expected outcomes. For large programs there exist software verification 
programs to facilitate the verification process. For a commercial solver it may not be possible to 
completely verify the program if the source code is unavailable. 


1.6.2 Validation 


Validation is really a check on whether the PDE is a good model for the real problem being studied. 
Validation means comparing numerical results with results from similar physical problems. Physical 
results may come from measurements from real life or from small-scale laboratory experiments. Either 
way, due to measurement errors, scaling problems and the inevitable failure of the PDEs to capture all the 
underlying physics, agreement between numerical and physical results will not be perfect and the user will 
have to decide what is ‘close enough’. 
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Introduction 


1.7 Exercise 1 


1. Assuming that t is time and x and y are spatial variables give the dimensions of the PDEs in a) to e) 


of Section 1.4. 


2. Classify the following PDEs: 


a) U, =2U,,, b) U,, +U,, =0, 


<< 


ak 


| “Ss 
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2. Fundamentals 


The finite difference method (FDM) works by replacing the region over which the independent variables 
in the PDE are defined by a finite grid (also called a mesh) of points at which the dependent variable is 
approximated. The partial derivatives in the PDE at each grid point are approximated from neighbouring 
values by using Taylor’s theorem. 


2.1 Taylor’s Theorem 


Let U(x) have n continuous derivatives over the interval (a, b). Then for a < x, x,th < b, 


Ua (x,) 


(n—D! +O(h"), (2.1) 


UGK +h)=U(K, HHU, (x, 4h? eh 


where, 


e U,(X,) is the derivative of U with respect to x evaluated at x = Xo. 


e Och") is an unknown error term defined in Appendix A. 


The usual interpretation of Taylor’s theorem says that if we know the value of U and the values of its 
derivatives at point x, then we can write down the equation (2.1) for its value at the (nearby) point x,th. 
This expression contains an unknown quantity which is written in as O(h") and pronounced ‘order h to the 
n’. If we discard the term O(h") in (2.1) (i.e. truncate the right hand side of (2.1)) we get an approximation 
to U(x,th). The error in this approximation is O(h"). 


2.2 Taylor’s Theorem Applied to the Finite Difference Method (FDM) 


In the FDM we know the U values at the grid points and we want to replace partial derivatives in the PDE 
we are solving by approximations at these grid points. We do this by interpreting (2.1) in another way. In 
the FDM both x, and x,th are grid points and U(x,) and U(x,+h) are known. This allows us to rearrange 
equation (2.1) to get so-called Finite Difference (FD) approximations to derivatives which have O(h") 
errors. Appendix A explains the meaning of O(h") notation. 
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2.3 Simple Finite Difference Approximation to a Derivative 
Truncating (2.1) after the first derivative term gives, 

U(x, +h)=U(x, #hU,(x,)+O(h?) (2.2) 
Rearranging (2.2) gives, 


U(x, +h)-U(x,) , O(h’) 
h h 


SERENE) 0(0) (byA.3.3) 


U,&)= 


Neglecting the O(h) term gives, 


U(x, th)-U(x,) 


U,(x,) ® 7 


(2.3) 


(2.3) is called a first order FD approximation to U ,(x,) since the approximation error = O(h) which 


depends on the first power of h. This approximation is called a forward FD approximation since we start 
at x, and step forwards to the point x,th. h is called the step size (h > 0). 


2.4 Example: Simple Finite Difference Approximations to a Derivative 


This simple example shows that our forward difference approximation works and has the stated order of 
accuracy. We choose a simple function for U. Let U(x) = x”. We will find the first order forward FD 


approximation to U , (3) using step size h = 0.1. From (2.3) the general first order forward FD 


approximation formula is, 


U(x, +h)-U(x,) 
h 


U,(x,) * (2.4) 


Substituting for U gives, 
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Replacing x, by 3 and h by 0.1 gives, 


2 aa 
U.Q) a ei 


The exact answer from basic Calculus is clearly U, (3) = 6 so the error in the approximation is 
6.1 —6=0.1. Repeating the problem with h = 0.05 (i.e. half the step size) gives, 


(3+0.05)?-3? _ 
0.05 


U,(3) 6.05 


The error is 6.05 — 6 = 0.05. The approximation formula (2.4) is first order so the errors should be 
proportional to h which is seen to be the case: halving the step size results in a halving of the error. 
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2.5 Constructing a Finite Difference Toolkit 


We now construct common FD approximations to common partial derivatives. For simplicity we suppose 
that U is a function of only two variables, t and x. We will approximate the partial derivatives of U with 
respect to x. As t is held constant U is effectively a function of the single variable x so we can use Taylor’s 
formula (2.1) where the ordinary derivative terms are now partial derivatives and the arguments are (t, x) 
instead of x. Finally we will replace the step size h by Ax (to indicate a change in x) so that (2.1) becomes, 


A 2 A n-l 
U(t, x, + Ax) = U(t, x,) + AU, (t, x,) 5 Ue (t,x,)+..4 Ug (to) + O(Ax") 


(n-1) 
(2.5a) 


Truncating (2.5a) to O(Ax’) gives, 
U(t, x, +Ax)=U(t, x,) + AxU,(t, x,) + O(Ax”) — (2.5b) 
Now we derive some FD approximations to partial derivatives. Rearranging (2.5b) gives, 


_ U(t,x,+Ax) -U(t, x,) O(Ax’) 

7 Ax Ax 

_ U(t, x,+Ax) -U(t, x,) 
Ax 


U,(t, x,) 


= UL(Ls,) 


x 


—~O(Ax)  (2.6a) 


Equation (2.6a) holds at any point (t, x,). In numerical schemes for solving PDEs we are restricted to a 
grid of discrete x values, x), X2, ... , Xn, and discrete t levels 0 = to, t), .... We will assume a constant grid 
spacing, Ax, in x, so that xj; = x; + Ax. Evaluating Equation (2.6a) for a point, (t,, x;), on the grid gives, 


U(t,, Xia) =U(tis X;) 


— O(A : 
AS O(Ax) = (2.6b) 


U,(t,, X;) = 


We will use the common subscript/superscript notation, 


U; =U(t,,x;) — (2.6c) 
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so that dropping the O(Ax) error term, (2.6b) becomes, 


U(x) 2-60 
Ax 
(2.6d) is the first order forward difference approximation to U, (t,,,X;) that we derived previously in 


approximation (2.4). We now derive another FD approximation to U,(t,,, x;). Replacing Ax by —Ax in 
(2.5b) gives, 


U(t, x, —Ax) = U(t, x,)-AxU, (t,x) +O(Ax”) (27a) 


Evaluating (2.7a) at (t,, Xi) and rearranging as previously gives, 


U?-U ae (2.7b) 


(2.7b) is the first order backward difference approximation to U,(t,,X;). 


Our first two FD approximations are first order in x but we can increase the order (and so make the 
approximation more accurate) by taking more terms in the Taylor series as follows. Truncating (2.5a) to 
O(Ax’), then replacing Ax by -Ax and subtracting this new expression from (2.5a) and evaluating at (t,, xi) 
gives, after some algebra, 


Una -U;?, (2.8) 


U,(t 
2Ax 


X;)* 


n? 


(2.8) is called the second order central difference FD approximation to U, (t,, x;) - 


We could construct even higher order FD approximations to U, by taking even more terms in the Taylor 


series but we will stop at second order approximations to first order derivatives. 


Many PDEs of interest contain second order (and higher) partial derivatives so we need to derive 
approximations to them. We will restrict our attention to second order unmixed partial derivatives i.e. U,, . 


Truncating (2.5a) to O(Ax*) gives, 
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2 3 


A A 
U(t, x, + Ax) =U(t, x,)+ AxU, (t, ie, Cat t=, alt x,)+O(Ax*) — (2.9a) 
Replacing Ax by -Ax in (2.9a) gives, 


Ax? 


“ay Uww (ts Xo) + O(AX") (2.9b) 


U(t, x, —Ax)=U(t, x,) — AxU, (t, x,) + “v,, (t,x,)- 
Adding (2.9a) and (2.9b) gives, 
U(t, x, + Ax)+U(t, x, — Ax)=2U(t, x, )+Ax?U,. (t,x,)+ O(Ax*) — (2.10a) 
Evaluating (2.10a) at (t,, x;) and using our discrete notation gives, 
Ut, +U? ,=2U"+Ax7U,, (t 


i+] i- 


x,)+O(Ax*) — (2.10b) 


n? 
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Rearranging (2.10b) and dropping the O(Ax’) error term gives, 


(2.11) is the second order symmetric difference FD approximation to U,, (t,, X;). These results are put 


into Table 2.1 to form a FD approximation toolkit. FD approximations to partial derivatives with respect 


eq UU 


U. 
ee (t,, X,)® = Ax? 


(2.11) 


to t are derived in a similar manner and are included in Table 2.1. 


partial derivative finite difference approximation type order 
OU oe 
ax ~ x U,,-U; forward first in x 
Ax 
OU 
— =U, U?-U? irst i 
ox i i-l backward first in x 
AX 
ou _ U n n 
ax ~ x Ua UE central second in x 
2Ax 
OU ror. U n n n 
a Oe Uj4-2U; +07, symmetric second in x 
Ax? 
ou — n+l n 
aet*~S Ui" -U; forward first in t 
At 
ou = n n-l 
aest~S U,-U; backward first in t 
At 
ou = UTt_un 
Ot ‘ cel ar central second int 
2At 
o-U n+l n n-l 
Ow Un Us" -2U;+U; symmetric second in t 
At? 


Table 2.1 Finite Difference Toolkit for Partial Derivatives 
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The above FD toolkit can be used to create a finite difference scheme (FDS) to obtain the approximate solution 
of a large number of PDEs simply by replacing each partial derivative by an appropriate FD approximation. 


2.6 Simple Example of a Finite Difference Scheme 


We construct a simple FDS to find an approximate solution of a simple PDE. This PDE will be studied in 
more detail in Chapter 4. For now it suffices to generate a simple FDS to provide motivation for further 
study. The 1D linear advection equation is, 


U,+vU,=0, (2.122) 


where the independent variables are t (time) and x (space). x is restricted to the finite interval [p, q] which 
is called the computational domain. v is a constant and the dependent variable, U = U(t, x). In addition to 
the PDE, we need initial conditions for U. Let the initial conditions be, 


UO, x) =f), p<x<q.  (2.12b) 


i.e. the initial value of U is given for every x value in the computational domain by a known function f(x). 


A solution to (2.12a, 2.12b) is a function U = U(t, x) which satisfies the PDE (2.12a) at all points x in the 
computational domain and all times t and the initial conditions (2.12b). U(t, x), the exact solution of 
(2.12a,b), is defined at an infinite number of values of the independent variables t and x. We will create a 
FDS to approximate U at a finite set of values of the independent variables. The approximate values of U 
on this finite set will be denoted by u. We proceed in stages. 


2.6.1 Step 1: Spatial Discretisation 


The computational domain (Figure 2.1) contains an infinite number of x values so first we must replace 
them by a finite set. This process is called spatial discretisation. 


Figure 2.1: 1D computational domain. 


For simplicity the computational domain is replaced by a grid of N equally spaced grid points. Starting with 
the first grid point at x = p and ending with the last grid point at x = q, the constant grid spacing, Ax, is, 


Ax=ld=P) 


oo (2.13a) 
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The values of x in the discretised computational domain are indexed by subscripts to give, 
Xi =p.%=ptAx,...,x=pt+(-l) Ax,...,xy=p+(N-1) Ax=q.  (2.13b) 
Since the grid spacing is constant, 
Xin =X: + Ax  (2.13c) 


The discretised computational domain is shown in Figure 2.2: 


Ax 
+ 
saaweies eae een eet agg Ge etcetera neg neta Ca age 
pExi X2 X3 Xn-1 Xn=q 


Figure 2.2 Discretised computational domain. 
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Fixing t at t= t, we approximate the spatial partial derivative, U, , in (2.12a) at each point (t,, x;) using 


the forward difference formula from the toolkit in Table 2.1 to give, 


n n 


Uy(ty xe 2.14) 
Ax 


Replacing U, in (2.1a) by its approximation (2.14) gives, 


U?, -U 
U,+v i+] i 
Ax 


0 (2.15) 
(2.15) is said to be in semi-discrete form since only the spatial derivative has been discretised. 


Note: The grid is also called the mesh and the operation of discretising the computational domain is called 
gridding or meshing. 


2.6.2 Step 2: Time Discretisation 


Fixing x at xX = x; we approximate the temporal partial derivative, U,, in (2.12a) at each point (tn, xj) using 


the first order forward difference formula from the toolkit in Table 2.1 (where At is the spacing between 
time levels) to give, 


U, & ur ur (2.16) 
On substituting (2.16) for U;, 2.15 becomes, 
ah ogee =0 (2.17a) 
At Ax 
which rearranges to give, 
UU *(u;,-U1) (2.17b) 
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Equation (2.17b) is an example of a FDS to approximate the solution of the PDE (2.12a). (2.17b) is a so- 
called time-marching scheme which enables U values at time level n+1 to be approximated from U values 
at the previous time level n. Since all U values are only known exactly at the initial time level (2.17b) is 
rewritten as, 


where u; is an approximation to U;=U(t,,x;,). 
Notes: 


1. (2.18) holds for each grid point x;, i= 1, 2, ... , N. u(ty, X;) is a numerical approximation to 
U(t,, Xi), the exact solution of the PDE (2.1a). 


2. u(0, x;) = U(0, x) but this will not be true in general for later times. 


3. uvvalues on the right hand side of (2.18) are all at time t, whereas on the left hand side u values are all 
at the next timed level t, + At = tn+ 


4. (2.18) is an example of a time-marching scheme in that (known) data for each grid point at time t, is 
used to find data at each grid point at the future time t, + At. This is called an iteration of the scheme. 
After an iteration of the scheme all u values at each grid point are known at time t, + At. These new 
values can be used as known data for another iteration of the scheme to give data for each grid point 
at the next time level. This process can be repeated until the required future time is attained. Iterations 
of (2.18) are performed by a computer program. 


5. The errors in approximating the spatial and temporal derivatives which are used to get (2.18) are 
O(Ax) and O(At) respectively and so (2.18) is said to be (formally) first order in space (x) and first 
order in time (t). 


6. The grid spacing, Ax, was determined by choosing the number of grid points, N. A larger N gives a 
smaller Ax and a (hopefully) more accurate solution as spatial derivatives are more accurately 
approximated. However as N increases compute time increases so there is a trade off between 
accuracy and speed. 
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7. The time step, At, is for the moment, chosen arbitrarily. However a smaller time step will mean that 
more iterations of (2.18) are needed to reach a stated future time which will obviously increase the 
compute time. In addition, since the result of each iteration is an approximation to the required 
solution, more iterations could cause the build up of more error. We will see later (Chapter 4 and 
Appendix C) that there is often a limit to the maximum size of a time step. 


8. (2.18) is said to be an explicit method since the value of u at the next time level is given by an explicit 
formula for each grid point. 


9. Later we will see that (2.18) doesn’t work for v > 0! There is more to FDS than meets the eye! 


2.7 Pen and Paper Calculation (very important) 


In practice numerical schemes are implemented by writing then running a computer program. Before 
doing this it is extremely useful to work through a pen and paper calculation for two reasons: 


1. To check understanding of the scheme. 


2. To be able to check results from the computer program against pen and paper results (verification). 
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Since we are doing a pen and paper calculation we will only use a small number of grid points. In a 
computer implementation of the calculation we would use many more grid points (hundreds, thousands, 
perhaps millions) for accuracy. We now set up a simple problem. 


Let p=0 and q=100, v=0.5 and let the initial conditions be, 


~0.01(x—45)? woe 
ui0.%{ ge VERS IO: 5) 


0, elsewhere 


Note that there is nothing special or realistic about these initial conditions. 


Let the (small) number of grid points be N = 11. Then by (2.13a), 


(100-0) _ 
ity: 


Ax 10 


The subscripts for the grid values go from | to 11 and are entered into the first row of Table 2.2. The actual x 
values of the corresponding grid points are 0, 10, 20, 30, ... , 90, 100 and are entered into the second row of 
Table 2.2. It remains to choose the time step At. Quite arbitrarily let At = 3. Then (2.18) becomes, 


u™=u"-0.15(u", ai?) (2.20) 


(2.20) is a FDS for calculating the solution to our problem at the next time level using data at the current 
time level. We start at time t,=0, i.e. n = 0, hence (2.20) becomes, 


u!=u!—0.15(u?,-u?) (2.21) 


Time level zero corresponds to the initial conditions. The initial u values are needed at the computational 


grid points. i.e. we need to know u) for i=1, 2, ..., 11. In general u” is only an approximation to the 


exact solution U(t,, x;) but at time to=0, u? = U(0, xj). 
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Using the initial conditions we get, 


-0.01(x,-45)2 
@ VIG)" 20<x< 70, 


0, elsewhere 


a =U(0; “| (2.22) 


Evaluating (2.22) at a few grid points gives, 
u°=U(0, x,)=U(0,0)=0, u2=U(0, x,)=U(0, 40)=e ">" =0.7788, etc. 


These initial values, u° , are entered into the third row of Table 2.2. 
1 


i | 1 2 3 4 5 6 7 8 9 | 10] 11 12 
x | 0 10 20 30 40 50 60 70 | 80 | 90 | 100 | 170 
ur | 0 0 0.0019 0.1054 0.7788 0 | 0 0 0 
u;' | 0 | -0.00029 -0.01363 0 
u2 


Table 2.2 Implementation of Finite difference Scheme (2.20) 
Having set up the initial data we use (2.21) to find u;' i.e. the u values at the next time level at each grid 
point. These new values are entered into the fourth row of Table 2.2. The first few values are found from 
the following: 
Putting i= 1 into (2.21) gives: 
1_.0 00 
ul=u?—0.15(u9—u' ) =0-0.15 (0-0)=0 


Putting i = 2 into (2.21) gives: 


ub=u9—0.15(u? —u5 )- 0 —0.15 (0.0019 — 0) = - 0.00029 
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Putting 1 = 3 into (2.21) gives: 
ub=u?—0.15(u? —u} )- 0.0019 — 0.15 (0.1054 - 0.0019) = -0.01363 


etc. 


However there is a problem! We are using forward differences for the spatial derivatives, i.e. to 
approximate the spatial derivative at a grid point we need the function value at the next grid point. This is 
OK for interior grid points but when we come to the last point on the right hand boundary of the 
computational domain (point index i= 11 corresponding to q = 100) we need data at point index i = 12 
which we DO NOT HAVE! In this case we have to invent a fictitious point with an associated function 
value. These points are called ghost points and the function values, ghost values. How we invent ghost 
values is based on the boundary conditions that define the particular problem we are solving. This topic is 
discussed in more detail in Appendix B. In our case we will assume that ghost point is at x,. = 110 and that 
the u value at this point takes the same value as its value at the nearest interior point (i.e. x:,) at all times. 


ie. U,=U,,.0=0,1,2,... (2.23) 
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Now we can use our numerical scheme to calculate uj, . 
Putting i= 11 into (2.21) gives: 


u!, =u? —0.15(u?, -u°,) = 0- 0.15 0-0) =0 


This completes the first iteration of the FDS (2.21) and row 4 of Table 2.2 and we have found the 
approximate solution at each grid point at 


t=t=At=3. 


Once all the values of u, have been calculated the same procedure can be used to find u; by a second 


iteration: the time indices in equation (2.21) are increased by | to give, 
2 1 1 1 
Uj; =U; ~0.15(u!,, —u; ) (2.24) 


and we repeat the previous process to fill in row 5 of Table 2.2 which is the approximate solution at each 
grid point at t= t= 2At=6. 


By successive iteration of (2.20), the solution can be found at each grid point at future time levels. At each 


iteration we use the known data at a particular time level to obtain the unknown data at the next time level. Of 
course this iterative procedure can be automated. A computer program for this scheme is given on the website. 


2.8 Exercise 2a 


1. U(x) =x’. Find the first order backward difference approximations to U, (3) using: a) h=0.1, 
b) h=0.05, c) h= 0.025 . 


2. Repeat Qla), b), c) using the central FD approximation and show that it is second order accurate. 
3. Following the text, derive the central FD approximation to a first order spatial derivative. 
4. Following the text, derive the symmetric FD approximation to a second order spatial derivative. 


5. Using exactly similar working for the spatial approximations derive all the time derivative 
approximation in Table 2.1. 


6. Complete Table 2.2 by carrying out pen and paper calculations (tedious but very important for 
checking your numerical algorithm). 


Download free ebooks at bookboon.com 


32 


Please click the advert 


Introductory Finite Difference Methods for PDEs Fundamentals 


2.9 Exercise 2b 


1. From the website download a computer program to implement (2.20) for the given problem. 


2. Read each line of the program and make sure you understand it. This program will be used as the 
basis for other programs later on. 


3. Verify the program by comparison to Table 2.2. 
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3. Elliptic Equations 


3.1 Introduction 


Elliptic PDEs form a class of PDEs that may be used to model steady state problems (i.e. the dependent 
variable remains constant over time). Solutions of elliptic PDEs are over closed regions on which boundary 
values are given in some way. These boundary values determine the solution of the PDE in the interior of the 
region. The two most widely used elliptic PDEs are Laplace’s equation and Poisson’s equation. 


In 2D, Laplace’s equation is: 


U,.+U,=0  .1) 


Laplace’s equation may be used to model a wide range of phenomena including steady state groundwater 
flow and temperature distribution over a region. Additionally Laplace’s equation can describe ‘potential 
flow’ which can be used in a simplified description of water flow amongst other things. 


In 2D, Poisson’s equation is: 


U,,+U,=f(x y). (3.2) 


Poisson’s equation may also be used to model a wide range of phenomena including gravitational fields, 
stress patterns and simplified viscous flow. 


The above PDEs can only be solved analytically for simple situations so we need to use numerical 
methods to obtain approximate solutions for cases of practical interest. In the following we focus on 
Laplace’s equation since it is simpler than Poisson’s equation and the techniques carry over easily. For 
simplicity the computational domain will be rectangular. 


3.2 Finite Difference Method for Laplace’s Equation 


The computational domain is discretised using constant grid spacings of Ax and Ay in the x and y 
directions respectively. Grid points are indexed by (i, j) in the usual way and the approximate value of U 
at grid point (i, j) is denoted by u,; . Figure 3.1 shows a rectangular grid with M and N grid points in the x 
and y directions respectively. u is known (=U) at the boundary grid points. It is required to find u at the 
interior grid points. 
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Figure 3.1: Computational grid showing interior grid points (black) and boundary grid points (white) 


Each partial derivative in Equation (3.1) is replaced by a symmetric FD approximation from our tool kit 
(Table 2.1) to give, 


1 


Wij 720; Uj + Us 720; FU = (3.3a) 
2 
Ax? Ay 


Letting b = Ax/Ay, (3.3a) can be rewritten to give, 


2 2 
a +b Uy jitb Uj ja (3.3b) 


a 2(1+b’) 


Uj j FU j 


Equation (3.3b) shows that u,; depends on its 4 surrounding values. This is called a 5-point stencil. 
Sometimes ‘compass notation’ is used and (3.3b) becomes, 


re u,tUy+b7u,+b7u, (3.3c) 
: 2(1+b’) 


where o denotes the current grid point and subscripts N, S, E and W denote its north, south, east and west 
neighbours respectively. 
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Notes: 


1. 


Using the indexing system in Figure 3.1 the unknown value of u nearest to the bottom left hand 
corner of the computational domain is uy. and the unknown value nearest to the top right hand corner 
of the computational domain is Uy.1.N-1. 


In an MXN grid there will be (M-2)x(N-2) unknown interior values of uj; which may be a very large 
number. 


Assuming that the boundary values of u are known then (3.3b) gives a system of (M-2)x(N-2) linear 
equations for uj; in (M-2)x(N-2) unknowns. 
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3.3 Setting up the Equations 


There are two basic methods of solving for u;;. Both methods set up a system of linear equations as 
follows. Letting ¢ = 1/(2(1+b’)), 


d=b*/(2(1+b’)), rearranging Equation (3.3b) gives, 


+du. 


i,j+l 


+cu, (3.4) 


i+,j 


u, j=CU;_, + du, 


Evaluating Equation (3.4) at successive grid points starting at 2,2 and sweeping along the rows first gives, 


Wy5 =CU, » +du,, + du; +CU,, 

U; =CU,, +du5, +du,, +CU 4 

Uys =Upso +dyy;, do Miz Toys (3.5) 
U3 =cu,, +du,, +du, 4 +CU;, 

Uvana =Uy-onat du M-I,N-2 + Uy in +CUM NA 


The boundary values uj 2, U2, Um2, Wi3, etc. are Known. In each equation the known u values are moved to 
the left hand side and the unknown u values moved to the right hand side (and positioned to preserve the 
ordering). For example in the first equation uj. and uy,; are known (being boundary values) so the 
equation becomes, 


—cu, ,—du, =— Un34 du, ,+cu, 5 


The result is that Equations (3.5) can be written as a single matrix equation, 
d=Au_ (3.6a) 

where, 

dis an (M-2)(N-2) by 1 matrix of known constants, 


A is a (M-2)(N-2) by (M-2)(N-2) matrix of known coefficients and 


u is a (M-2)(N-2) by | matrix of unknowns and 
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u- (U2 2 U3,.2--- Uaw-1),2 U2,3 U3 3... UdW-1),3 -- - Ueno) 
The solution to (3.6a) may be written symbolically as, 
u=A'd_ (3.6b) 


As A may be very large we must study efficient ways of finding u. 


3.4 Grid Convergence 


Before looking at solution methods we make the following important point which applies to all grid-based 
numerical schemes. Clearly the accuracy of the numerical results depends on the size of the computational 
grid. Ideally we would like a grid converged solution 1.e. a solution that does not change significantly 
when more grid points are used. Grid converged solutions can be studied formally by grid convergence 
indices but for us it is enough to compare numerical solutions with more and more grid points until there 

is no significant difference. A good strategy is to implement a scheme using a small number of grid points 
(for verification of the program) then successively double the number of grid points until numerical results 
don’t change perceptibly. In this way we can be confident that inaccuracies in the numerical solution are 
not caused by the grid. We present other questions about the accuracy of schemes in Appendix C that is 
best read later. 


3.5 Direct Solution Method 


One way to achieve a solution of (3.6a) is by using standard Gaussian elimination. This is a so-called 
direct method. 


Note: 


We state again that in our convention the grid is traversed from left to right starting at the bottom, i.e. we 
start at position (2, 2) go along the whole row, repeat for the next row etc. and eventually end at position 
(M-1, N-1). 


The direct method is illustrated by an example on a small grid (much too small for a real computation but 
useful for verification of a program). Consider the following 5 by 4 grid of u values where 
Ax = Ay, 
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Figure 3.2 Boundary (white) and interior (shaded) u values on a 5x4 grid. 
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The u values at the boundary grid points are known, e.g. u;,; = 6.1. We want to find the unknown values of 
u at the interior grid points, i.e. Uz, Us,2, U42, 23, U33, U43. Since Ax = Ay, Equation (3.4) simplifies to, 


+u,; 


i,j+l 


+U; 


i+L,j 


u,=(u,.,40 \/4 (3.7a) 


ij ijl 
In compass notation this is written, 
u,=(uy tus tuy+u, \4 (3.7b) 


Starting on the second row, evaluation of (3.7a) at grid point (2, 2) gives, 


U22 = (Ui. + U2,1 3 U23 + U32 y/4 = (7.2 + 6.8 + U23 + U3.2 y/4 (3.8a) 


Evaluation of (3.7a) at grid point (3, 2) gives, 


U3.2 = (U22 rU31 7 U3 3 U42 \/4 = (U2 2 + 7.74 U33 Uu42 y/4 (3.8b) 


Evaluation at grid point (4, 2) gives, 


U42 = (U3,2 + U4 + U43 + Us2 /4 = (U3 + 8.7 +43 + 9.4/4 (3.80) 


Moving to the next (third) row, evaluation of (3.7a) at grid point (2, 3) gives, 


U23 = (ui3 U2,2 U24 U3 3 \4 = (3.4 U22 + 8.9 4 U3,3)/4 (3.8d) 


Evaluation of (3.7a) at grid point (3, 3) gives, 


U33 = (U2,3 + U3.2 + U3 4 + U43 y/4 = (U2 3 + U3.2 + 8.9 + u43)/4 (3.8e) 


Finally evaluation of (3.7a) at grid point (4, 3) gives (fill in this for yourself), (3.8f) 


These equations can be written out in standard form as 6 simultaneous linear equations in which the order 
of the unknown variables is as above i.e. U2 2, U3,2, U4, U23, U33, U4 3. 
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Rewriting Equation (3.8a) in standard form gives, 
-14/4 =-1 Wo. + (1/4) 32 + 0 gn + (1/4) 13 + 033 + 0 U3 
Rewriting Equation (3.8b) in standard form gives, 
-7.7/4 = (1/4) ung -1 us. + (1/4) uy. + 0 un3 + (1/4) 33 + 0 43 


Equations (3.8c-f) are rewritten similarly. The 6 simultaneous linear equations are then written as the 
matrix equation, 
—14/4 -1 1/4 0 1/4 0 OY 422 
—7.7/4}) |1/4 -1 1/4 0 1/4 O | Ys2 
PoP ob ob GG du 
= meh 4a oy 
U3 


U3 


U4 3 


which can be solved by standard Gaussian elimination. For large systems this isn’t very efficient so we 
will study efficient iterative methods next. 


3.6 Exercise 3a 

1. Finish labelling all the grid points in Figure 3.1. 

2. Check that (3.3a) is correct by referring back to the FD toolkit in Table 2.1. 
3. Derive (3.3b) from (3.3a). 


4. Given that Ax = 4 and Ay = 2, write down the expression for u3,4 using both subscript and compass 
notation. 


5. Write down Equation (3.8f) in the space provided in the notes. 
6. a) Write down Equations (3.8c-f) in the standard way. 
b) Hence complete the matrix equation (3.9). 


7. Solve (3.9) by Gaussian elimination and check your solution by back substitution. 
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3.7 Iterative Solution Methods 


Equation (3.9) can be expressed as, 


Au=b_ (3.10) 
For practical problems A is likely to be a large matrix which makes the direct solution of (3.10) 
computationally inefficient. More efficient methods use iterative approaches where an initial estimate for 
u 1s updated to form a better estimate. The process is repeated until the distance between successive 
estimates is less than some pre-defined tolerance (assuming that the iterative process converges to the 
solution). Since the output from each iteration is a vector of u values we define what is meant by distance 
between vectors as follows. 
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3.7.1 Distance between Vectors 


Let x =(x), X2, ... , Xn) and y =(y}, yo, ... , yn) be two vectors in R™. The distance between x and y is 
denoted by d(x , y) and defined by, 


d(x,y)=max(|x;-y; |)=-y|,  G.1)) 
(this is sometimes called the ‘infinity norm’). 


In the following we illustrate three variants of the iterative approach with respect to the test problem of 
Figure 3.2 where Ax = Ay. 


When Ax = Ay, our 5-point formula is, 


+u. 


i,j+] 


+u, 


i+],j 


w= (Uti V4 G.12) 


3.8 Jacobi Iteration 


We introduce the iteration index as a superscript, m, and write (3.12) as the Jacobi formula (also called the 
point-Jacobi formula), 


m+l_{..m m m m 
u =(u", tu" tu", +u",, yl 4 (3.13) 


(Note that this formula assumes Ax = Ay). 


For each interior grid point (i, j), ujj at the next iteration (m+1) is found from (3.13). Once an iteration has 
been completed for all interior grid points we compute the distance between vectors u™”! and u™. If, 


m-+1 m | 


lu = —u |,<tol, (3.14) 


+ 


where tol is a pre-defined tolerance, the iterations terminate and the solution to (3.10) is u" ' otherwise the 


iterations continue. 
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To start the iteration (at index 0) values must be given for the unknown interior values of u;;. These values 
can be set to zero or interpolated from the known boundary values. The following example of Jacobi 
iteration applies to the previous ‘test problem’ with all interior starting values set to zero and a tolerance 
of 0.5 x 10°. 


3.8.1 First Iteration 


Equation (3.13) with m=0 gives, 


1 _(.0 
ul =(u? tu 


0 


0 
a +u, 


0 
+u,; i+1,j 


i,j+l 


ya 


Evaluating at each interior grid point gives, 


uh =(u?,+u8 +09 ,+u2, )/4=(7.2+6.8+0+0)/4=3.500 
uh =(u8 +0? +09 +9, )/4=(0+7.7+0+0)/4=1.925 , similarly 
0 =4.525, 0 4=4.325, 0, 4=2.225.0,,=4.525.. 6.15) 


Having completed the first iteration we test for convergence. i.e. we measure how close w’ is to u’. 
u’ = (0, 0, 0, 0, 0, 0), u’ = (3.5, 1.925, 4.525, 4.325, 2.225, 4.525), 


therefore, 
{u'—u"|,,=max (| 3.5— 0 |,|1.925—0],| 4.525— 0 |,|4.325 — 0 |,|2.225 — 0|,|4.525 —0]) 


This is greater than the specified tolerance of 0.5 x 10° so the iteration is repeated to find wu’ etc. 
Eventually (we hope!) after p iterations, 


lu’—u? “|, <0.5x10° 


and the iterative procedure terminates. The solution is u’. 


A program to implement Jacobi iteration is available from the website. 
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3.9 Gauss-Seidel Iteration 


This is a potentially more efficient version of Jacobi iteration. We note that in Equation (3.13) some of the 
updated u;; values are already available for use in the iteration formula. In our way of traversing the grid 
when we reach grid position (i, }) we have already updated u;.,; and uj;.; therefore we can use these values 
in the 5-point formula which becomes, 


m+l__ m+l1 m+l m m 
Le =(uteuttleum tut, (4 (3.16) 


(Note that this formula assumes that Ax = Ay). 
This is called the Gauss-Seidel formula (also called point-Gauss-Seidel) and the implementation is the 


same as for Jacobi. As an example of Gauss-Seidel iteration we repeat the previous problem using starting 
u values of zero (and where we don’t yet have an updated u value we use its current value). 
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3.9.1 First Iteration 


Equation (3.16) with m=0 gives, 


a tue, 4 


1 (1 1 
Uu. =u -+U. +Uj i iti 


i,j i-l,j i,j-l 


Evaluating at each interior grid point gives, 


iia (ul, +u,,+u5,+U55 ) 4. 
Note that we don’t have Usa yet so we use Uy 5501, 100 give, 
=(u?, +02, +u9,+u9, )/4=(7.24+6.8+0+0)/4=3.500 
1,2 2,1 233 3,2 " 7 * . 
is (u), +u,,+U;,+U4> y/ 4 
Note that we now have ee but not Wea SO we use Wl is to give, 
=(u) +0? ,+u9 tu, )/4 =(3.5+7.7+0+0)/4=2.8 
2,2 3,1 3,3 4,2 7 * 
a 4= (ul, +4, +U4,+Us, y/ 4 
Note that we now have Wi6 but not U4 1,80 we use U4)» to give, 


=(u!,+u9+u) +02, )/4 =(2.8+8.740+9.4)/4=5.225 
(3.17) 


The remaining values are calculated similarly. 
Notes: 
1. We are assuming that the iterative procedures converge — this needs further study. 


2. For simplicity we have taken Ax = Ay. It is a simple matter to generalise the iterative formulae for the 
case of different grid spacing in the x and y directions (start from Equation (3.4)). 
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3. 


It is quite tricky to write a general program for an M by N mesh. The difficult part is to map the 
(M-2)x(N-2) array of interior grid values to a system of (M-2)(N-2) linear equations in (M-2)(N-2) 
unknowns. 


3.10 Exercise 3b 


6b. 


Check the direct solution to (3.9) by using Scilab (or Matlab). 
Find |x—y |, for x=(1,3,5,8,4,-2,4,0,2), y=(2,-4,4,1,4,2,-3,0,2). 


Using pen and paper perform a second iteration for each unknown for the Jacobi test problem 
(equations (3.15)). 


Write a program to implement Jacobi iteration for the test problem above and verify by comparison 
to your pen and paper calculations (and by comparison to the direct solution). Comparison programs 
can be found from the website. 


. Using pen and paper complete the first iteration of the Gauss-Seidel calculations for the test problem. 


Perform a second iteration. 


. Write a program to implement Gauss-Seidel iteration for the test problem using zero starting values 


and verify by comparison to 5a,b. 


For a given tolerance compare the number of iterations taken using Gauss-Seidel and Jacobi methods 
for the test problem. 


Write a Jacobi iterative solver for a general rectangular grid and validate it on the test problem 
(a program is available on the website for comparison). 


Repeat Q7. using a Gauss-Seidel iterative solver (a program is available on the website for 
comparison). 


3.11 Successive Over Relaxation (SoR) Method 


The idea behind this method is that in an iterative formula the point value at the new iteration depends on 


the old point value plus some error (residual) at that point. 


Le = up =uFRS 6.18) 
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Ri is the difference between successive iterates of uj; and is called the residual. It might be possible to 


speed up the convergence of the iterative scheme by weighting the residual on the right hand side of (3.18) 


appropriately. We write, 


uy, =ur+wRF, (3.19) 


w is called a relaxation parameter. For 0 < w < 1 (3.19) corresponds to under-relaxation and for 
1<w<2 (3.19) corresponds to over-relaxation. This idea is applied to improve the point-Gauss-Seidel 
method. The point-Gauss-Seidel iteration formula (for Ax = Ay) is,’ 


m+l__/{,.m+l m+l m m 
i= (oa a y 4 (3.20) 


which can be re-written as, 


m+l1 m+l m 
-ijtUija 


i,j+] 


Te =u) +(u —4u7 HO tO \ 4  (3.20a) 
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which is the of the same form as (3.18) with, 


m m+l m+l m 
Rj j= =(u i- 4,j+Uj 140; Faaltee 


+u™ )/4 (3.206) 
Hence we may convert this to, 


u™=(1—w)u™+w(uramiun tum )/4 — B.21) 


Ui; i,j-l ijt Wij 


which is called the (point) SoR method. This method is implemented in the same way as the previous 

iterative methods. If w = 1 (3.21) reduces to the Gauss-Seidel method. The question is what is the best 
choice of w for fastest convergence? This is a difficult question to answer in general and we must use 

numerical experiments to find an approximate best value. For the interested reader some convergence 

analysis for the three iterative methods is given in Appendix D. 


3.12 Line SoR 


As we shall see next, the use of a classical tridiagonal matrix method can greatly improve the efficiency of 
the previous SoR method. This is because we can update the solution along a whole line of grid points at 
once instead of simply point by point, hence the name /ine SoR. The line of grid points will normally be a 
whole row or column in the grid. We start from the finite difference toolkit approximation to Laplace’s 
equation in compass notation (3.3c) 


m+1 


m+] _ u, r+uy +b?(uy+ug ) 
7 2(1+b’) 


(3.22) 


where we have introduced the iteration index m and assumed the usual ordering working through the grid 
points left to right and bottom to top. (3.22) is the Gauss-Seidel method (3.16) in compass notation on a 
rectangular mesh with b = Ax/Ay. The corresponding point SoR method is, 


ws" (wou +59 farsa stages] G23) 


which is the compass notation form of (3.21) on a rectangular mesh. We proceed by moving the East and 
West terms onto the left hand side, with the East term written at the m+1 level, i.e. (3.23) becomes, 


201+ b? ue" — wut" —wut"=b, (3.24) 
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Where 
b, =2(1+b’)U—w)ut + wb?(ut +ug") (3.25) 


m+1 


and noting that the data u, on the right hand side is known if we assume a bottom to top ordering for 


the calculation. (3.24) is now in tridiagonal form. The corresponding matrix equations for the unknown 
data u at the interior grid points (i = 2, M-1) ina single row j are: 


(20462) —w Pee | |. Bg 
-w 2(1+b*) -w 0 Tea by: 
=| - (3.26) 
0 —w 2(1+b’) -w ' , 
[ —w = 2(1+b*) funy} | bx, 


where we see that a tridiagonal matrix is one in which the entries in the main diagonal and in the diagonals 
above and below it are in general non-zero, with the remaining entries zero. A tridiagonal system like 
(3.26) can be solved very efficiently using the Thomas algorithm (see website) for the m+1 iterated values 
of u along one complete line or row of the mesh at once. This is the /ine SoR method, or SoR by lines. For 
convenience, the matrix equations (3.26) assume zero end-point boundary values. Dirichlet or von 
Neumann boundary conditions (see Appendix B) involve only marginal changes to the entries in the 
tridiagonal matrix and the right hand side vector b. (3.26) also assumes that the nodes in each row j are 
numbered left to right from i=1 to M where at the end points (i=1 and i=M) boundary conditions are 
applied. Instead of visiting each node in the mesh at each iteration, we solve a complete row at a time by 
solving (3.26), working up the rows bottom to top to complete one iteration cycle. The procedure is 
illustrated in Figure 3.3. The row sweeps are continued until the solution values converge to the required 
accuracy. In order to maintain diagonal dominance of the equations (3.26), and retain computational 
efficiency, we must ensure that w < 1+b’(on a square mesh w<S 2). 
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U=Up=const 


U=Up=const 


| sweep by rows 
bottom to top 
Apply Thomas 
algorithm then 
advance to next 
row up 


U=U9=const 


U=Up=const 


Figure 3.3 SoR by Lines 


The faster convergence of line SoR compared to standard (point) SoR is due to the greater influence of the 
boundary values that affect all nodes at each sweep. It can be shown that the optimum value of the 
relaxation parameter w is given by the smaller root of, 


t’w’-l6w+16=0 (3.27) 


where t = cos(x /(M- 1))+ cos(x I(N- 1)) and M and N are the number of grid points in the x and y 


directions respectively. 
3.13 Exercise 3c 


1. Write a general program to implement SoR and validate it on the test problem. Experiment with 
tolerances and choice of relaxation parameter. 


2. Repeat Q1 using line SoR. 
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4. Hyperbolic Equations 


4.1 Introduction 


Hyperbolic equations describe many time-dependent (transient) phenomena (e.g. fluid flow) and are 
characterised by a speed of propagation of information (often called ‘wave speed’). Consequently future 
solution values can only be affected by past values in a limited neighbourhood. Hyperbolic equations may 
also admit discontinuities in the solution variables requiring advanced numerical treatment. A proper 
mathematical treatment of hyperbolic equations involves the study of characteristics that is beyond the 
scope of this book. Instead we use experimental results to illustrate key concepts. Attention is focused on 
the 1D linear advection equation that, although not strictly hyperbolic according to the classical definition, 
contains the essential elements of hyperbolic PDEs and, as we will see later, is a component of the shallow 
water equations that are truly hyperbolic. We will look at numerical results from the simple FD scheme in 
Chapter 2 to solve the 1D linear advection equation and this will lead to a series of basic questions about 
numerical schemes in general. 
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4.2 1D Linear Advection Equation 


The linear advection equation (also known as the transport equation) may be used in a model of various 
phenomena like the movement of pollutant in a river or the movement of an air/water interface inside a 
sophisticated numerical method (e.g. Volume of Fluid method). We look at the simplest version of this 
equation. In one spatial dimension the linear advection equation is: 


U,+vU,=0  (4.1a) 


Here the independent variables are t (time) and x (space). The dependent variable U is a function of t 
and x. 


A PDE on its own does not give a complete model of a particular problem. In addition to the PDE, we 
need to give the correct initial and/or boundary conditions for the dependent variable for the particular 
problem. Let the initial condition for U in (4.1a) be, 


U(0, x) = f(x) (4.1b) 


i.e. the initial value of U is given over the spatial domain by a known function f(x). For the moment we 
will assume that the spatial domain is infinite so that we can ignore any boundary conditions. 


4.2.1 An Interpretation of the Linear Advection Equation 

If we interpret (4.1a, 4.1b) as a (partial) model of the transport of a soluble pollutant by a 1D river then 
U(t, x) is pollutant concentration at time t and position x along the river and v is the (constant) velocity of 
the river. Equation (4.1b) gives the initial pollutant concentration at each point along the river. For a fixed 
value of t the graph of U against x is called the concentration profile. Suitable units measure t in seconds, 


x in metres and U in kg/m. A dimensional analysis of (4.1a) shows that v is measured in m/s which is 
correct for velocity. 


A solution to (4.1a, 4.1b) is a function U = U(t, x) which satisfies (4.1a) and the initial conditions (4.1b) at 
all points (x, t). 


4.2.2 Exact Solution of the Linear Advection Equation 
It can be shown that the exact solution to (4.1a, 4.1b) is, 


U(t, x) =f(k-vt) (4.1c) 
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This means that U(t, x) is just the initial concentration profile, f(x), translated vt metres along the x axis. 
For v > 0, the translation is to the right and for v < 0, the translation is to the left. In either case the 
pollution moves downstream at the speed of the river. The following diagrams illustrate this using an 
arbitrary made up initial concentration profile. 


concentration U(0, x) = f(x) 
U (kg/m) 


2 4 x (metres) 


Figure 4.1a Initial concentration profile. 


concentration wo, 
U (kg/m) 


) 4 2+vt Atvt -x (metres) 


Figure 4.1b Concentration profile after time t (solid) compared to initial profile (dashed), v > 0. 


Figure 4.1a shows a triangular initial concentration profile along the river. Concentration is only non-zero 
in the interval 2 <x <4. At time t the concentration profile is the same shape as the initial concentration 
profile but has been translated vt units downstream. Physically equation (4.1a) models pollutant transport 
in the absence of diffusion so the pollutant is just carried along at the speed, v, of the river. This model is 
unrealistic but is useful for learning purposes. It is important to note that the model is NOT a model of the 
flow in the river which flows with constant velocity v (we will look at models of water flow later as they 
are more complicated). 


The linear advection equation is a good PDE to study because it has an exact solution and we can assess 
the performance of FD schemes by comparison to this. 
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4.3 Results for the Simple Linear Advection Scheme 


Graphical output from the simple FD scheme of Chapter 2 is given for a range of parameters. In all cases 
the same initial condition profile given by Equation (2.19) was used and the computational domain, 

[0, 100], was discretised by N = 101 equally spaced points. There is a brief discussion of the output in 
each case which will motivate the analysis later in the chapter. 


4.3.1 Test Case 1 


comparison of solutions to du/dt + v du/dx = 0, + numerical, o analytical 
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Figure 4.2a Comparison of numerical (+) and exact solutions (0) to the 1D linear advection equation using 
first order forward differences in both space and time using v = 0.5, At = 0.3, 10 time steps. 


The simulation was run for 0.3x10 = 3 seconds. The initial condition profile has moved about 1.5 metres 
to the right. This is what we would expect since v = 0.5m/s is positive and distance equals speed 
multiplied by time. There is ‘good’ agreement between numerical and exact solutions. 
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4.3.2 Test Case 2 


comparison of solutions to du/dt + v du/dx = 0, + numerical, o analytical 


1.2 T T T 


concentration u 


T T T T 


Figure 4.2b Comparison of numerical (+) and exact solutions (0) to the 1D linear advection equation using 
first order forward differences in both space and time using v = 0.5, At = 0.3, 25 time steps. 
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The simulation has been run for 0.3x25 = 7.5 seconds. The numerical concentration peak has moved to the 
right place but is higher than the exact solution. More worryingly there is some noticeable divergence 
from the exact solution in the numerical solution around x = 15 and x = 67. Something is going wrong! 


4.3.3 Test Case 3 


comparison of solutions to du/dt + v du/dx = 0, + numerical, o analytical 


2 0 T T T T T T T T T 
4 + 


15} +s | 


concentration u 
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Figure 4.2c Comparison of numerical (+) and exact solutions (0) to the 1D linear advection equation using 
first order forward differences in both space and time using v = 0.5, At = 0.3, 44 time steps. 


The simulation has been run for 0.3x45 = 13.5 seconds and numerical results have gone haywire. Note 
that the vertical scale has changed and the numerical results have ‘blown up’. Something is terribly wrong 
with this scheme! 
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4.3.4 Test Case 4 


comparison of solutions to du/dt + v du/dx = 0, + numerical, o analytical 
1 T ‘ ¥ | T T T T T 


concentration u 
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Figure 4.2d Comparison of numerical (+) and exact solutions (0) to the 1D linear advection equation using 
first order forward differences in both space and time using v = -0.5, At = 0.3, 44 time steps. 


The simulation has again been run for 13.5 seconds but the sign of the velocity, v, has been changed. As 
expected the concentration profile moves to the left. Results look reasonable this time so there is a lack of 
symmetry with Test Case 3. Obviously the behaviour of the numerical scheme is influenced by the sign 
of v. 
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4.3.5 Test Case 5 


comparison of solutions to du/dt + v du/dx = 0, + numerical, o analytical 
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Figure 4.2e Comparison of numerical (+) and exact solutions (0) to the 1D linear advection equation using 
first order forward differences in both space and time using v = -0.5, At = 13.5, 1 time step. 
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The simulation has again been run for 13.5 seconds with v = -0.5, as before. The time step has been 
increase to 13.5 so that the end result is achieved in a single iteration. As before the concentration profile 
moves to the left. However results look bad and there is evidence that the numerical solution has started to 
“blow up’! The time step, At, could be too big. 

A number of questions arise naturally from the results of these numerical experiments: 

Q1) Can we design other FD schemes to get more accurate results? 

Q2) How do we know the numerical results won’t ‘blow up’ at some future time? 


Q3) How do we know we have an accurate solution (in the absence of an analytical solution)? 


These questions are addressed by the theory in Appendix C which should be read now. 


4.4 Scheme Design 


We will design FD schemes to solve the 1D linear advection equation (4.1a). First we introduce some 
useful compact notation. 


4.4.1 Operator Notation 


FA) ral 2 
Let, 0, denote a’ O,, denote re O,,, denote 


Then (4.1a) can be written, 
0,U+vo0,U=0 (4.2a) 


.6,U=-vd,U — (4.2b) 


By definition, 


6,U=0,(0,U)  (4.3a) 
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Using (4.2b), (4.3a) gives, 
0,,U=0,(-vo VU) (4.3b) 


=—vo, (0. U) 


=—vo,(—vo, U) 
=v'd,.U (4.3c) 


Now we can design some FD schemes to solve Equation (4. 1a). 


For fixed x, the Taylor expansion of U(t + At, x) to order 3 gives, 


2 
U(t+ At, x)=U(t, x) +AU, +—U,+ 0(At) (4.4a) 


which, in operator notation is, 


2 
U(t+ At, x)=U(t, x)+Ata,U +—0,U+ O(At) (4.4b) 


Using (4.2b) and (4.3c) the partial derivatives with respect to t can be replaced by partial derivatives with 
respect to x to give, 


2 
U(t+ At, x)=U(t, x)-Atvd,U +2 -v90,,U+ O(At®) — (4.4c) 


Let, 


cy 
Le (At)=1—Atvo, ts? Y O (4.5) 


XX 
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Then (4.4c) can be written, 
U(t+At, x)=L, (At)U(t, x)+O(At’) (4.6) 


L, (At) is a differential marching operator. To design FD schemes to solve (4.8a) we simply redefine 
L, (At) by replacing each continuous partial derivative by a finite difference approximation (denoted by 
5, 5xx) to give, 


2 


L,(At)=1-atvs, + v5 (4.7) 


XX 


L, (At) is now a difference marching operator. Different FD approximation choices for 5,, 5, give rise to 


different FD time marching schemes. The general (second order in time) FD time marching scheme for the 
1D linear advection equation (4.8a) can now be written as, 


u™t=L (At)u" — (4.8a) 
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Written out in full this is, 


n+l n n At? 2 n 
u; =U; “EO ee u; — (4.8b) 


XX 1 


(As previously defined, u;’ is the approximation to the exact solution U(t, xi) at the i” grid point and the 


n"" time step). Some examples of FD schemes are now given. 


4.4.2 Example 1: Forward Time Centred Space (FTCS) Scheme 


n n 
: U,,,;—U,_ 
From our FD toolkit we choose 6,u; =—“1—" and 6, u; =0. 
x 


Equation (4.8b) becomes, 

n+l n vAt n n 

i FU Uj “Uj (4.9) 
7 ut) 

This is the FTCS scheme and it has the following stencil. 


time level 


n+1 


: ZN 


i-1 i i+] 
+ spatial steps 
2Ax 


Figure 4.3 Stencil for the FTCS Scheme 
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Notes: 


1) The scheme is first order in time and second order in space (see Appendix C for definition of the order 


of a scheme). 


2) Ghost values are required at both left and right ends of the computational domain (see Appendix B for 
boundary conditions). 


4.4.3 Example 2: First Order Upwind (FOU) Scheme 


n n 
u; —U;_; 


From our FD toolkit we choose 6, u; = and 6,,u; =0 . Equation (4.8b) becomes, 


X 


uur“ (up—up) (4.10) 


This is the FOU scheme and it has the following stencil. 


time level 


n+1 


7 


¢——> spatial steps 


Figure 4.4 Stencil for the FOU Scheme 


Notes: 
1) The scheme is first order in time and first order in space. 


2) A ghost value is required at the left end of the computational domain. 
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4.4.4 Example 3: Lax-Wendroff Scheme 


u 
From our FD toolkit choose 6, u;’ = 


Equation (4.8b) becomes, 


(4.11) 
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This is the Lax-Wendroff scheme and it has the following stencil. 
time level 


n+1 


: aa 


i-1 i it] 
spatial steps 


Figure 4.5 Stencil for the Lax-Wendroff Scheme 


Notes: 
1) The scheme is second order in time and second order in space. 


2) Ghost values are required at both left and right ends of the computational domain. 


4.4.5 Example 4: Lax-Friedrichs Scheme 


This is the same as the FTCS scheme except that the first term on the right of (4.9) is replaced by the 


: ; . ae UU; tu; ; 
average of its 2 neighbouring values, i.e. u; is replaced by eae . As in the FTCS scheme we choose 


n 


n 

U;,—U 

oS se and 6, u;=0. 
xX 


Equation (4.8b) becomes, 


n 


n 
ot Ui FU vAt (u’ 


i 5 AAx* —U; ) (4.12) 
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This is the Lax-Friedrichs scheme and it has the following stencil. 
time level 


n+1 


: a Ss 


i-1 i i+] spatial steps 


Figure 4.6 Stencil for the Lax-Friedrichs Scheme 


Notes: 
1) The scheme is first order in time and first order in space. 
2) Ghost values are required at both left and right ends of the computational domain. 


3) Although this scheme appears to be quite similar to the FTCS scheme its performance is very 
different (see Appendix C). 


4.5 Multi-Level Scheme Design 

So far all our schemes have been based on using data at the current time level (n) to advance to the next 
time level (n+1). This approach can be extended to multi-level schemes by performing Taylor 
approximations at t - At and using algebraic manipulation as we shall now see. Replacing At by — At in 


(4.4a) gives, 


A 2 
U(t—At, x)=U(E, x)-AtU, +5-U, +O(At’) (4.13a) 


Subtracting (4.13b) from (4.13a) and gives, 


U(t+ At, x)-U(t — At, x)=2AtU,+O(At?) — (4.13b) 
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In operator notation this is, 
U(t + At, x)-U(t — At, x)=2Atd,U+O(At*) — (4.13c) 


Using (4.2b) this becomes, 
U(t + At, x)-U(t — At, x)=- 2vAtd,U+O(At*) — (4.13d) 


Dropping the error term, replacing the differential operator by the difference operator and using the usual 
discrete notation gives the general FD scheme, 


u;=u?'—2vAté,u (4.14) 
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4.5.1 Example 5: Leap-Frog Scheme 


UU. : 
In (4.14) we choose 6, u; =—*—“" to give, 
2Ax 


n n- vAt n n 
Vie =u; 1 ie (u2,-u’, ) (4.15) 


This is the Leap-Frog scheme and it has the following stencil. 


spatial steps 


time level 
n+1 
At 
n 
n-1 
Figure 4.7 Stencil for the Leap-Frog Scheme 
Notes: 


1) The scheme is second order in time and second order in space. 


2) Ghost values are required at both left and right ends of the computational domain. 


3) Initial conditions are required at two time levels! 


4.6 Exercise 4a 


la. Download the simple linear advection solver ‘linearadvection’ (from Chapter 2). 


1b. Run ‘linearadvection’ on Test Cases 1-5 and check that the output agrees with the figures in the notes. 
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lc. Run ‘linearadvection’ with different run times, time steps, spatial steps and values of v (both positive 
and negative). In each case write down in your own words how the numerical scheme behaves in 
comparison to the exact solution. 


ld. Are there any values of At, Ax, v, ntimesteps that give good results? 


2a. Copy ‘linearadvection’ to ‘linearadvectionFOU’. Change the solver in this new file so that it 
implements the FOU scheme. Verify your code by comparison to a pen and paper calculation as per 
Table 


2.2. Important: the FOU scheme requires a left hand ghost point and hence in the first non-ghost point has 
index 2 and the last non-ghost point has index N+1. 


2b. Run ‘linearadvectionFOU’ using the parameters from Test Cases 


1-4. In each case write down in your own words how the numerical scheme behaves with reference to the 
analytical results. 


2c. Run ‘linearadvectionFOU’ with different run times, time steps, spatial steps and values of v (both 
positive and negative). In each case write down in your own words how the numerical scheme 
behaves with reference to the analytical results. 


2d. Are there any values of At, Ax, v, ntimesteps that give good results? 
2e. Show experimentally that the FOU scheme is first order in space. 


3. Copy ‘linearadvection’ to ‘linearadvectionLW’. Change the solver in this new file so that it 
implements the Lax-Wendroff scheme and verify the code. Repeat 2b-e for ‘linearadvectionLW’. 


4. Doall exercises in Appendix C and relate theoretical stability results to your experimentally derived 
results. 


4.7 Implicit Schemes 


The previous schemes are called explicit schemes because data at the next time level is obtained from an 
explicit formula involving data from previous time levels. This leads to a (stability) restriction on the 
maximum allowable time step, At (see Appendix C). Now we come to a different type of scheme — 
implicit, in which data from the next time level occurs on both sides of the difference scheme that 
necessitates solving a system of linear equations. There is no stability restriction on the maximum time 
step which may be much larger than an explicit scheme for the same problem. In implicit schemes the 
time step is chosen on the basis of accuracy considerations. 
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4.7.1 Example 6: Crank-Nicolson Scheme 


This is an implicit scheme. In previous examples spatial derivatives are approximated at time level n. 
However values change between time level n and time level n+1 so a better approximation to spatial 
derivatives could make use of data at both time levels. Let, 


n n n+l n+l 
Sits ing) 0<a<i: 416) 
2Ax 2Ax 


This is a weighted average of central difference approximations at times levels n and n+1. Choose 
6,,U; =0. Then for a = % (4.8b) becomes, 


n+l 


n n n+l 
yay? vies Ui 4 Bis Ui (4.17) 


' * QAx 2 2 
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This is the Crank-Nicolson scheme and it has the following stencil. 


time level 


n+1 


At 


i-1 i it] 


spatial steps 


Figure 4.8 Stencil for the Crank-Nicolson Scheme 


Notes: 

1) The scheme is first order in time and second order in space. 

2) Ghost values are required at both left and right ends of the computational domain. 

3) The scheme is implicit so values at time level n+1 are found by solving a system of linear equations. 


4.7.2 Implementation of the Crank-Nicolson Scheme 


At 
Letting c=~—-, (4.17) is, 
Ax 


n n n+l n+l 
yttay?— C} Ui UA ‘ Ui Ui (4.18) 
2 2 


Re-writing (4.18) gives, 


n+] __ n n n n+l n+l 
4u;"=4u; —c(u’, =1. 70, =U ) (4.19a) 


i+] i-l 
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Rearranging so that data from the same time level is on the same side gives, 


n+l 


—cu;, +4u 


n+l n+l 


n n n 
* +cu;,, =cu;,+4u;—cu;,,  (4.19b) 


The data at time level n is assumed known so (4.19b) is simplified by replacing the right hand side by d; 
to give, 


cur +4u™—cut!=d? — (4.19¢) 


i i+l 
For each grid point i= 1, 2, ... , N, we write out (4.19c) to give, 


n+l 


1 1 
—cu, +4u;" +cu;” =d; 
—cur+4u3t + cust! =d, 


n+l n+l n+l n 


—cu, +4u;° +cu, = (4.194) 


n+l 


n+l n+l _ jn 
—cuy.,+4uy,+cuy =di, 


n+l n+l n+l 


—cuy_,+4uy, +cuy,, =dy 


(4.19d) is a system of N linear equations in what looks like N+2 unknowns! However u?*! and ut! on 


the left hand side of (4.19d) are ghost values which may be known directly or can be calculated in terms of 
neighbouring values depending on the type of boundary condition given in the problem (see Appendix B). 
In the former case we can move these known values to the right hand side of (4.19d) and, letting 

dt =d"+cu?"',d\=d cut", (4.19d) becomes, 


n+l 


4u""'+cus =d? 
n+l n+l n+l _ yn 
=cu; +40," + cu; =d> 
—cus'+4uh!+cuy"! de 
(4.19e) 
n+l n+l n+l _ qn 
=O, 6 4u, cu =di1 
n+l n+l _ an’ 
—cuy,+4uy =d 


This system is expressed as the matrix equation, 


Au™'=d"  (4.19f) 
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where, 


4 0 
-c 4 c 0 
0 -c 4 

A=| 0 -—c 
0 0 
0... 


The solution to (4.19f) is obviously, 


nt] 


0 
0 yt a, 
0 .. 0 i a, 
ec 0... Of,u™’=| : | d"=] : 
eee dna 
—c¢ 4 Cc we lee 
0 -c 4 
=A'd" (4.19g) 
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Notes: 


1. (4.19g) is solved at each time step and the solution updated iteratively. 

2. Inpractical problems A may be very large (e.g. 1000 x 1000) so an efficient matrix inversion method 
may be needed (see Chapter 3). 

3. In this example boundary values are known so A is constant and needs only to be inverted once. For 
problems where the left hand side boundary values are calculated in terms of neighbouring values 
(e.g. Derivative boundary conditions, see Appendix B) the first and last rows of A will vary from 
time step to time step. 

4. Inour example A has a special structure — it is tridiagonal. This special structure permits the use of 
the efficient Thomas algorithm for matrix inversion algorithm (see downloadable code for Chapter 3). 


Results for the usual test case are given in Figure 4.9. 


Comparison of Crank-N numerical (+) and exact solutions to du/dt+vdu/dx=0 
1 T T T T T T T an T T 


0.8; 


0.4; 


concentration, u 


0.2+ 


-0. 2 L 1 L | L | 1 | i 
0 10 20 30 40 50 60 70 80 90 100 


distance, x 


Figure 4.9 Comparison of numerical (+) and exact solutions (0) to the 1D linear advection equation using 
the Crank-Nicolson scheme with v = 0.5, C = 2.0, 15 time steps. 
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4.8 Exercise 4b 


1. 


Use von Neumann stability analysis (Appendix C) to show that the Crank-Nicolson scheme is 
unconditionally stable. 


Using the Crank-Nicolson scheme for the 1D linear advection equation with p = 0, q = 100, 


v =0.5, N=5, C= 1 and Dirichlet boundary conditions uj =0=u,,, Vn, and the usual initial profile. 


Write down the system of linear equations for the first time step of the Crank-Nicolson scheme and 
express them as a matrix equation. Invert the (5x5) matrix (by pen and paper or use a package) and 
hence find the concentration values at each grid point after the first time step. 


Using your calculation from Q2 to verify your program and write a program to implement the Crank- 
Nicolson scheme using N = 50 grid points. Experiment with different time steps and compare your 
numerical results to the exact solution. Write a short report detailing the characteristics of the scheme 
(tip: use the code from a previously written scheme as a basis for your new code). 


What do you see? 


A skateboarder? 


We see an All Options Junior Trader leaving 
the office after a successful day. 
He has the freedom to make his own decisions. 
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5. Parabolic Equations: the Advection-Diffusion 
Equation 


5.1 Introduction 
The advection-diffusion equation belongs to the class of parabolic PDEs. In | spatial dimension it is, 
U,t+vU, =K,U,, (5.1) 


K, is called the diffusion coefficient (in the x direction). If K, = 0 then (5.1) is the linear advection 
equation which we studied in Chapter 4. Using our previous interpretation of the linear advection equation 
in which U = U(t, x) is river pollutant concentration and v is the speed of the flow, (5.1) is a more realistic 
description of pollutant transport. Not only does the initial pollutant move downstream with velocity v, the 
pollutant also diffuses into the surrounding water at rate K, (the presence of second order spatial 
derivatives often indicates a diffusive process). 


Figure 5.1 illustrates a pure advective process (K = 0 in (5.1a)) compared to an advection-diffusion 
process (K > 0 in (5.1b)). As we have seen previously for pure advection there is no change in the initial 
concentration profile as it moves downstream at speed v. When diffusion is introduced the initial 
concentration profile moves downstream at speed v but also spreads out (diffuses) and reduces in height 
over time. It is important to realise that Figure 5.1b shows the exact solution — the spreading and reduction 
in height is a consequence of the underlying physical process. Of course we have seen similar behaviour 
with numerical solutions hence we talk of ‘numerical diffusion’ that is a feature of most numerical 
schemes. 


pure advection: initial profile (dotted line) and 2 later solutions advection-diffusion: initial profile (dotted) and solution at later times, 
T r T r r r r r r T T T T T 


0.9 4 0.9 | 
0.8 al 0.8 7 
0.7 4 0.7 7 

3 I. 4 3 

5 0.6 é 0.6 J 

£05 4 Eos 4 

i= i— 

oO oO 

2 3 

8 0.4- 4 S$ 0.4 4 
0.3 7 0.3 4 
0.2 4 0.2 | 
0.1 4 04 | 

(0) 1 L L L L L L L L 0 y L + 1 L 1 1 1 
0 10 20 30 40 50 60 #70 80 90 100 0 10 20 30 40 #50 60 70 480 490 #100 
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Figure 5.1 Time evolution of solutions for advection (a) and advection-diffusion (b) 
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5.2 Pure Diffusion 


To keep things simple we will start with the pure diffusion equation, 


U,=K,U 


A ~U,  (.2a) 
5.2.1 Analytical Solution of the Pure Diffusion Equation 
Equation (5.2a) has an analytical solution for special conditions that is useful for checking numerical 


schemes. It can be shown that over the x interval [0, 1] with initial condition, U(0, x) = sin(zx) and 
boundary conditions, U(t, 0) = U(t, 1) = 0 for all t, the solution of (5.2a) is, 


U(t,x)=e *"'sin(ax) — (5.2b) 


5.2.2 Finite Difference Scheme for the Pure Diffusion Equation 


We start simply and use our FD toolbox to replace the time derivative by a first order forward difference 
and the spatial derivative by a symmetric difference so that (5.2a) becomes, 


n n n 
U, Yin 2U; +0; 
At a‘ Ax? 


(5.3a) 


This can be rearranged to give the explicit FD scheme, 


AtK 
upt=upeo (u2,—2u’4u",) (5.3b) 


n 
i+] 


Note that the presence of u;,, and u,;', indicates the need for left and right ghost values. Equation (5.3b) 


is written compactly as, 


ul =(I- 2r)u?+runtiu", — (5.3¢) 


- 7 i+] 
where, 


AtK , 
t= ; 
Ax 


(5.4) 
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Since this scheme is explicit the time step, At, will be limited by stability constraints. A von Neumann 
stability analysis (Appendix C) gives, after some manipulation, 


G=(1-—21r)+2rcos(k Ax) (5.5a) 


where G is the amplification factor. For stability, |G |<] which gives, 


2 


Atc= (5,56) 
2K, 


The problem with the analytical solution (5.2b) was run for 300 time steps with 51 grid points using 
scheme (5.3b) with K, = 1. The maximum time step given by (5.5b) was multiplied by a safety factor of 
F=0.9. Figure 5.2 shows the results. It can be seen that the initial profile reduces in height but does not 
change its location. This is what we expect from a pure diffusion problem. The numerical and exact 
solutions are close. 


pure diffusion: initial profile (--) and numerical (+) and exact solutions 


0.9 


0.8 


0.7 


conggatration u 


0.5 ye Xt 
0.4 
0.3 
0.2 


0.1 


Figure 5.2 Comparison of exact and numerical solutions for a 
special pure diffusion problem. 
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5.2.3 Exercise 5a 
1. State the units of K, given standard units for the other variables in (5.1). 


2. Given the initial concentration in Figure 5.1 draw plausible solutions to (5.1) at two later times on the 
same graph when v = 0. 


3. Show that (5.2b) satisfies the special initial and boundary conditions. 

4. Show by partial differentiation that (5.2b) is a solution of (5.2a). 

5. Obtain (5.3b) from (5.3a). 

6. Obtain (5.3c) from (5.3b). 

7. By conducting a von Neumann stability analysis obtain (5.5a) then (5.5b). 


8. Write a program to solve (5.2a) using the scheme (5.3c). Check your program on the special diffusion 
case and reproduce Figure 5.2. 
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5.3 Advection-Diffusion Equation 


We discretise the advection-diffusion equation (5.1) using a first order forward difference for the time 
derivative, a first order backward difference for the first space derivative (assuming v >0) and a (second 
order) symmetric difference for the second space derivative to give, 


n+l n n n n n n 
u, —u, U. =U We 4H 20. Hu, 
i i +V 1 i-l =K, i+l i-l (5.6a) 
At Ax Ax 
which can be rewritten as, 

n+l n vAt n n K At n n n 
—|u;'—u;), }+—-lu 3, — 2u;) tu; 5.6b 
1 Ax ( 1 1 :) Ax? ( +1 :) ( ) 


(5.6b) is an explicit FD scheme for solving the advection-diffusion equation. It is first order in time and 
space. As usual we need to find the allowable time step by a stability analysis. For convenience (5.6b) is 
written, 


ult=Su" +(-R-2S)u"+(R+S)u", — (5.6c) 


where, 


R=——. (5.7a) 


K 
S=—=—. (5.76 
Ag (5.7b) 


By von Neumann stability analysis it can be shown that, 
1-R-2S20 (5.8) 
from which it follows that, 


Ax? 


M<—_— 43) 
vAx+2K, 
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Note that when v = 0, (5.9) reduces to the previously calculated stability limit for pure diffusion (5.5b) and 
when K, = 0 (5.9) gives the well known stability limit for pure advection. Results from scheme (5.6b) are 
given in Figure 5.3 for v = 0.5, K, = 0.1, N = 101 at time t = 57s where At is its maximum value in (5.9) 
multiplied by a safety factor of 0.9. The simulation required 45 time steps. By inspection the peak of the 
concentration profile has moved about 28m. This correct as v = 0.5 and 0.5 x 57 = 28.5m. The peak value 
of the profile has decreased and its width has increased as would be expected from a diffusive process. 


advection-diffusion: initial profile (dotted) and unsplit numerical solution (+) at later time, 
1 T T t T T t T T T 


0.9 - a 


0.8 - 


0.7 - 


0.6 - 


0.5 - 


0.4+ 


concentration u 


0.3 - 


0.2; 


Figure 5.3 Initial profile and numerical solution (+) to the 1D advection-diffusion equation at t = 57s. 
Important: 


When running a simulation always check your results as much as you can. There are various ways to do 
this: 


Check 1. If there is an exact solution to a special problem run the simulation on it and compare numerical 
and exact results. 


Check 2. Do a pen and paper calculation on a small problem and compare results with your numerical 
output. 


Check 3. Subject your results to a critical scrutiny — are they what you expect to get given the physical 
process you are modelling? 
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Of course the whole purpose of a simulation is to model some aspect of reality so the only real check (of 
the model) is to compare it against real data. 


5.4 Exercise 5b 


1. Derive (5.6a) 
2. Obtain (5.6b) from (5.6a). 
3. Obtain (5.6c) from (5.6b). 


4. Starting with (5.6c) carry out a von Neumann stability analysis on the solution and obtain (5.8) and 
hence (5.9). 


5. Calculate Ax and hence At for the simulation whose results are shown in Figure 5.3. 


6. Draw an approximate sketch of Figure 5.3 if the simulation had been run for 70 seconds. 


t : t 
A —o ,b) lim —2 
Ax>0 At. 


7. For v >0 show that, a) iim =|. What do you conclude from this? 
x 


8. Write a program to implement scheme (5.6b) and reproduce Figure 5.2. (hint: extend code from a 
previous scheme). 
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6. Extension to Multi-dimensions and Operator 
Splitting 

6.1 Introduction 

So far we have looked at PDEs in one (spatial) dimension. The real world is 3D so we must extend schemes to 


multi-dimensions (it should be noted however that there are many useful 1D and 2D models of 3D processes). 
For simplicity and to illustrate key concepts we will consider the 2D linear advection equation, 


U,+v,U,+v,U,=0  (6.1a) 


Like the previous 1D linear advection equation this is another so-called ‘model’ equation in that it is a 
simplified version of reality. As for the 1D linear advection equation we can interpret (6.1a) as a (partial) 
model of pollutant transport in a 2D river. U= U(t, x, y) is the pollutant concentration at time t at position 
(x, y) in the river (plan view) and v, and v, are the (constant) components of river velocity in the x and y 
directions respectively. Given initial conditions, 


U(O, x, y) = g(x y), 
it can be shown that (6.1a) has the exact solution, 
U(t, x, y) =g(k-vx tp y—vyt) (6.1b) 


i.e. the concentration profile is simply translated (advected) along at the speed of the river without 
changing shape (there is no diffusion term in (6.1a) — it models pure advection). 


6.2 2D Scheme Design (unsplit) 
Our first approach to 2D scheme design is to directly discretise the PDE. These schemes are said to be 


unsplit as opposed to split that we will look at later. We will design FD schemes to solve (6.1a) in an 
exactly similar way to Chapter 4. In operator notation for partial derivatives (6.1a) can be written, 


0,Ut+v,0,U+v,0,U=0  (6.2a) 


“.0,U=-v,0,U-v,0,U — (6.2b) 
By definition, 


0,U=0,(0,U) (6.2c) 
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Using (6.2b, 6.2c) gives, 
0,,U=0,(-v,0,U-v,0,U) (6.2d) 


=—v,0,(0,U)-v,0,(0,U) 
=—v,0,(0,U)—v,0,(@,U) 


=-v,0,(-v,0,U-v,0,U)-v,0,(—v,0,U-v,0,U) 
=V, 0. ,U+Vv,v,0,,U+v,v,0,,U+v, “Oy yU 


xy yx y xX 7 Xy 


=V, "0. ,U+2v,v,0,,U+v, "On U_ (6.2e) 


Bey AY: 


Now we can design some FD schemes to solve the 2D linear advection equation. Holding x and y constant, 
the Taylor series for U(t+At, x, y) expanded about U(t, x, y) to O(At’) is, 


2 


At 
U(t + At, x,y)=U(t, xy)tAtd, U+——8,U +O(At*?) — (6.3a) 
Using (6.2b, 6.2e) to replace the time derivatives in (6.3a) and rearranging gives, 


+2v,v,0,,+V20,, JU +O(At*) — (6.3b) 


xX XX xX y Xy y yy 


U(t + At, x, y)=U(t, x,y)-At(v, 0, +v,0 putea 
Let, 


At? 
Lyy (At)=I-At(v,0, +v,0 dei ae t2V,VyOqy +V50y,) 


(6.4) 
Then (6.3b) can be written, 
U(t+At, x,y)=L,, (At) U(t, x, y)+ O(At*?) — (6.5) 


Lyy (At) is a differential marching operator and is second order in time. To design FD schemes to solve 
(6.1a) we simply redefine L,,, (At) by replacing each continuous partial derivative by a FD approximation 


(denoted by 6,, dy etc.) to give, 
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2 
Lyy (At)=Il-At(v,6, +, 8, (v28 +2V,.V/OytV,Oyy) (6.6) 


XXX y~ xy 


This is now a 2D difference marching operator. Different FD choices for 6,, 5,, etc. give rise to different 
FD time marching schemes. 


The computational domain is 2D so we need a grid of points in x and y directions. We assume a rectangular 
computational domain with constant grid spacing of Ax and Ay in the x and y directions respectively. At is the 
time step. Extending our usual notation by using subscript j for the index in the y direction let, 


uj = Ulta, Xi, yj) (6.7) 


The general 2D FD time marching scheme for the 2D linear advection equation (6.1a) can now be 
written as, 


uj; =Lyy(At)u;, (6.8a) 
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Written out in full this is, 


ij X ~ XX xX" y~ xy 


At? 
u;; =u;,-At(v,6, +5, JUL + (v8 +2v,V,6,y+V,5,,)u;; — (6.8b) 


Some examples of FD schemes are now given. 


6.2.1 Example 1: First Order Upwind (FOU2D) Scheme 


n 


n n n 
u;.—u?,. a yn, 
j 9 -l, . jo 
From our FD toolkit we choose, 5, u;,=—2— 1 Sul,= iar a 
3 


x dsj 


6,, Uj ;=6,,U;,=5,,U; =0 . Equation (6.8b) becomes, 


ce antes xy ij 


ue <0 (ah 8,.)=6 (a ey ) (6.9) 


ij i,j 


where, 
v, At v,At _, . 
C= i C,= i are the Courant numbers in the x and y directions respectively. This scheme is 
x y 


simply the FOU scheme applied to each dimension. 
Notes: 
1. The scheme is first order in time and first order in each spatial dimension. 


2. Ghost values are required at the left side and the lower side of the computational domain. This is 
shown in Figure 6.1 for a grid with 5 points in the x direction and 3 points in the y direction. 


3. The 2D FD scheme (6.9) could simply have been obtained by direct replacement of the partial 
derivatives in the PDE (6.1a) using approximations from the FD toolkit. Our scheme is an example of 
an unsplit scheme because it is obtained directly from the 2D PDE. 


4. The scheme uses backward differences for the spatial derivatives and is therefore only an upwind 
scheme for positive velocities. When one or more velocities are negative the scheme will fail unless 
forward differences are used appropriately. 
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5. The scheme is explicit and will be subject to a stability constraint on the allowable time step. In the 
1D case the heuristic stability approach reasoned that information could not travel across 2 grid 


points in a single time step giving the condition, At<Ax i\yv| . In 2D it seems reasonable to take, 


At<Fmin (ax / 


Vx 


,Ay/lv,|) (6.10a) 


where F < 1, is a ‘safety factor’ which can be determined from numerical experiments if necessary. In 
principle it is possible to undertake a stability analysis for a 2D scheme in the same way as for the 
corresponding 1D scheme. This is algebraically complicated. A von Neumann analysis gives, 


C.+C,<1  (6.10b) 


which is a very restrictive condition on At. 


i —_> «<> 
Ax 


Figure 6.1 2D Computational mesh showing real grid points (black) and ghost grid points (white) 
for the FOU2D scheme 


Code to implement the FOU2D scheme (6.9) may be downloaded from the website. Graphical output from 
the above scheme is given in Figure 6.2. 
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exact solution to du/dt + vx du/dx + w du/dy = 0 contour plot of exact solution 


concentration 


0 20 40 60 80 100 


FOU2D scheme solution to du/dt + vx du/dx + w du/dy = 0 


concentration 


Figure 6.2 Comparison of exact and numerical solutions for the (unsplit) FOU2D scheme 
6.2.2 Example 2: Lax-Friedrichs Scheme in 2D 


From our FD toolkit we choose, 


n n n n 
U.,.—U:,. Uu...—U:. 
n _ ~i+l,j i-l,j n _ ~i,j+l i,j-l n n n 
OS pee he Big ae 
n n n n 
Uy FY ia FU ja FU 


and estimate u;', by (i.e. the mean of surrounding values). Equation (6.8b) 


4 


becomes, 


n 


n n n 
nt Yin Ui FU HU C, (u® age i Cy (u" = ) (6.11) 
uj = 4 3 Ming Bi) 5 iin 7 Fi . 
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Notes: 
1. The scheme is first order in t and in both x and y. 
2. Ghost values are required at all sides of the computational domain. 


3. This scheme is less diffusive than the FOU2D scheme and also works for positive and negative 


velocities. 


4. A von Neumann stability analysis gives, 


Cy+Cjs1/2 (6.12) 


which is a very restrictive time step condition. 
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6.2.3 Example 3: Crank-Nicolson Scheme in 2D 


By an obvious extension of the 1D scheme (6.8b) becomes, 


n n n+l n+l n n n+l n+l 
uty? C, (a Ui) j 4 Ui j es) _ a Ui i a Ui iat =) (6.13) 


2 2 

Notes: 

1. The scheme is first order in t and second order in both x and y. 

2. Ghost values are required at all sides of the computational domain. 

3. The scheme is implicit so values at time level n+1 are found by solving a system of linear equations. 

4. Von Neumann stability analysis shows that this scheme is unconditionally stable. 

5. Implementation of the Crank-Nicolson scheme involves solving MN linear equations in MN 
unknowns at teach time step where MN is the number of computational grid points. The matrix can 
be written as a penta-diagonal matrix. For large MN the computational effort at each time step may 
be prohibitive unless an efficient matrix inversion algorithm is used. 


6.2.4 Conclusions 


From the above analysis we can see that FD schemes extend naturally from 1D to 2D (and 3D). The main 
issues are: 


1. Itis often difficult to perform stability analysis. 
2. For explicit schemes the allowable time step may be prohibitively small. 
3. For implicit schemes the computational effort at each time step may be prohibitively large. 


4. It would be nice to be able to exploit the relative simplicity and larger time steps of 1D schemes for 
2D problems. The following analysis shows that we can actually do this and more! 
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6.3 Operator Splitting (Approximate Factorisation) 

This method replaces a single scheme to solve a complicated PDE by a sequence of simpler schemes 
which solve related PDEs and which together solve the original PDE (up to a specified order of accuracy). 
Operator splitting can be used to split up a PDE by dimension or by components or a combination of both. 
We begin by looking at dimensional splitting. 

6.3.1 Dimensional Splitting 

The basis of this method is to split our N dimensional PDE into N, 1D PDEs and use a sequence of 1D FD 


schemes to solve the problem. For simplicity we will look at the 2D linear advection equation (6.1a). We 
split the (6.1a) into the following 2 PDEs (which are both 1D in space), 


U,+v,U,=0  (6.14a) 
U,+v,U,=0  (6.14b) 


From our 1D work in Chapter 4 we know that a (second order in time) differential marching operator for 
solving (6.14a) is, 


ae 
De ee ae O (6.15a) 


X~ XX 


Similarly a differential marching operator for solving (6.14b) is, 
AG 2g 
Ly (At)=l-Atv,0,+ Vy yy (6.15b) 


Consider the operator sequence, 
Ly(At) Lx(At) (6.16) 
We show that this sequence is ‘the same’ as Lxy defined in (6.4) which we know solves the 2D PDE 


(6.1a). By ‘the same’ it is enough to show that the difference is of an appropriate order of At. Multiplying 
out (6.16) gives, 


A APs 
Ly (ADL x(At)=(I-Atv, 0,+——v¥0,, )0-Atv,0,+-VyOx.) (6.174) 
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2 
=l—Ai(v,¢, +V,0,)4+ (20,429 v,0O +v\0,,)+O(At*) (6.17b) 


xy RY 


=Lyy(At)+O(At*®) — (6.17c) 


Hence the sequence of 1D operators, Ly(At) Lx(At), is the same as the 2D operator Lxy(At) up to O(At?). 
This shows that the split scheme, 


u=L,(At)[Lx(Atju?,] (6.18) 
can be used to solve (6.1a) instead of the unsplit scheme (6.8a). It is now just a matter of replacing the 


differential operators, by appropriate difference operators like we did before: e.g. Lx becomes the 
difference marching operator, 


2 
Ly (At)-I-Atv, 8+ 


°8 


X7Xx * 


6007 'A'G Sweysds Wy] 282] @ 


It’s only an 
opportunity if 
you act on it 
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ev 
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Notes: 

1. (6.18) has the same temporal order as (6.8a). 

2. The time step in (6.18) is the minimum allowable time step for each 1D operator and may be much 
larger than the time step for the unsplit 2D scheme (6.8a) which is one of the advantages of 
dimensional splitting. 

3. (6.18) is implemented by first applying the Lx operator to each row of the grid. This produces new 
values of the dependent variable at each grid point. The Ly operator is then applied to the new data 
along each column of the grid. 

4. The split operator sequence (6.18) is one of many split operator sequences which solve (6. la). 


5. Lx and Ly may be totally different FD schemes (although they should be the same order). 


6. Our approach is called ‘dimensional splitting’ because a 2D PDE (6.1a) has been split into 2, 1D 
PDEs (6.14a, 6.14b). 


6.3.1.1 Example 4: Split FOU2D Scheme 


The FOU scheme is used for both Lx and Ly difference operators. 1.e. (for positive velocities) we choose, 


n n n n 
5 ul= ii Jini 5 u2 Mah Mahe 
u 90, U; i= 


x4i = , 0,,U;,=6,,u; =0 . So that, 
, Ax 


xx Yj yy “i,j 


L, (At)=I-Atv,8,, Ly(At)=I-Atv,6, . 


Neglecting terms of O(At’) the split sequence (6.16) is, 
Ly (At)L, (AD=(I-Atv, 8, )(I-Atv,8,) 
=LAt(v,6,+v,6, ) 


Equation (6.18) can be written, 


X sweep: uij=u",—C,(u",-u",,] (6.19a) 


y sweep: i =Uij =6 fi, -095) (6.19b) 
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where, 


Ui, j (‘u bar’) is the new data at the grid point i,j produced by the Ly operator and C, = re C= A 
x y 
are the usual Courant numbers. 


Notes: 


1. (6.19a) corresponds to uij=Lx(At)u?, and (6.19b) corresponds to OF SL (Ai)ing. 


; x A Ax A . : 
2. At<min cases where —— 2 are the maximum time steps for the Lx and Ly operators 
vl Wy v 
y 


xX xX 


; 


respectively. These expressions were found previously by von Neumann stability analysis of the 1D 
FOU scheme. 


Code to implement the split FOU2D scheme (6.19a, 6.19b) may be downloaded from the website. 
6.3.1.2 Operator Splitting for Stiff Problems 


The allowable time step, At, for an explicit scheme depends on the ‘flow’ velocity and the grid spacing in 
each direction. Ideally we want to time-march with as large a time step as possible for computational 
speed. This cannot always be achieved by the ‘standard’ splitting (6.18) as we will see from the following 
example. Let the allowable time steps in x and y directions be At, and At, respectively. For the split 


A 
FOU2D scheme, At=min(At, At, ) where At,=—~ : 
Vv 


xX 


A 

At =") Now suppose that |v,| >> |v,| (>> means 
| 

“much greater than’) then, assuming that the grid spacings are of the same order, 


At, << Aty. 


This defines a ‘stiff problem: the time scales for each component of the problem are very different. The 
‘standard’ split operator sequence, 


L,(At) Ly(At) (6.20) 
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requires that each operator marches using the same time step which has to be the minimum of At, and At, 


for stability. Clearly the time step for Ly is being constrained by the much smaller allowable time step in 


the x direction. It might be more efficient if we could change the splitting so that Ly could march in steps 


of At,. We can! 


Suppose that nAt, = At, where n is an integer. i.e. At, is n times as small as Aty. Consider the new split 


operator sequence, 


Ly (At, [Ly (At, /n)J" (6.21) 


In this sequence Ly is marching at its maximum time step, At,, and Lx is marching at (1/n) of At, so to 


‘catch up’ we need to use Lx n times. The following analysis shows that this approach works. 
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From the split FOU2D scheme example, 


L, (At)=I-Atv,6,, Ly (At)=l-Atv,6, . Substituting into (6.21) and using the Binomial expansion gives, 
(I-At, v6, )(I-(At ,/n)v,6,)” 
=(I-At, v,3, )(-n(At, /n)v,3, +0(At?)) 


(6.22) 
=(I-At, v,8, )(I-At, v,8, +O(At?)) 


2 
=I-At,(v,8, +v,8,)+O(At?) 
which is Ly(At,) Lx(At,) up to O(At’). This shows that the sequence (6.21) is valid. 
Notes: 


1. IfAt, is not exactly a whole number times as small as At, replace At, by nAt, where n is the largest 
integer such that nt, < At,. 


2. (6.21) says apply the operator Lx to each row of the grid n times using a time step At,/n each time and 
using the updated values of u at each application. Then apply the Ly operator to each column of the 


grid once using a time step At,. 
3. Since the multiplication in (6.22) is commutative another valid operator sequence is, 
[L, (At, /n)]/ "Ly (At, [L, (At, /n)]”. 
4. IfLyand Ly are second order (they aren’t in our example) and n is even then the symmetric sequence, 
[L, (At, /n)]"?Ly (At, [Lx (At, /n)]"” 
is also second order. Where possible it is almost always better to use symmetric sequences. 
6.3.1.3 Exercise 6a 
1. Show that (6.1b) is a solution to (6.1a). 


2. Go through the calculations to get (6.2e). 


3. Derive (6.3b). 
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4. Write out (6.8a) in full to get (6.8b). 
5a. Download ‘linearadvectionFOU2D’ and make sure you understand every line of code. 
5b. Run ‘linearadvectionFOU2D’ and describe the output for the given parameters. 


5c. By gradually increasing the safety factor F and running the code find its maximum value for a stable 
scheme based on the heuristic time step formula. 


5d. Find the maximum time step for the FOU2D scheme based on the von Neumann stability analysis 
result and alter your program to use this time step. 


5e. Change the sign of the x component of velocity and run the code. What happens to the numerical 
solution and why? 


i) Fix the problem so that your code can deal with v, < 0 and v, > 0. Test your program. 
ii) Fix the problem so that your code can deal with any velocities. Test your program. 


6. Find the expression for the time step based on von Neumann stability analysis for the 2D Lax- 
Friedrichs scheme to solve the 2D linear advection equation. 


7. Copy ‘linearadvectionFOU2D’ to ‘linearadvectionLF2D’ and change the solver in this new file so 
that it uses the 2D Lax-Friedrichs scheme (don’t forget ghost values) and repeat QS. 


8. By multiplying out and collecting terms, derive 6.17c from 6.17a. 


9. Show that Lx(At) Ly(At), is the same as Lxy(At) up to O(At’) and hence write down another operator 
sequence to solve the 2D linear advection equation. 


10. Show that for the 2D linear advection equation with equal velocities in x and y directions, the time 
step for the LxLy split FOU2D scheme is twice that for the unsplit FOU2D scheme when grid 
spacings in x and y directions are equal. 


11. Compare time steps in Q10 when the Lax-Friedrichs scheme is used for all operators. 


12. Using the notation of (6.19ab) write down the FD scheme to solve the 2D linear advection equation 
with the Lax-Friedrichs scheme for all operators using the LxLy sequence. 
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13. Copy a suitable previous code to ‘linearadvectionFOU2Dsplit’ and modify it to solve the 2D linear 
advection equation using the split operator sequence, Ly(At) Lx(At), using the FOU scheme for each 
operator. 


14. Copy a suitable previous code to ‘linearadvectionFOU2Dsplit’ and change it to solve the 2D linear 
advection equation using the split operator sequence [L, (At/2)]L, (At)[L, (At/2)] using the Lax- 


Friedrichs scheme for each operator. 
6.3.2 Term Splitting 
PDEs may contain several terms corresponding to different physical processes. As an example we use the 
1D advection-diffusion equation (5.1). Rather than solving the advection-diffusion equation directly it 
may be more efficient to use term splitting to solve each part separately and combine solutions. The 


following time step analysis shows why this could be a preferred option. 


Assuming v > 0, using a first order forward difference for the time derivative and a first order backward 
difference for the first space derivative, by previous analysis the time step for pure advection is, 


At,SAx/v_— (6.23a) 
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Using a first order forward difference for the time derivative and a symmetric difference for the space 
derivative, by previous analysis the time step for pure diffusion is, 


2 


Ax 
At, <— ___ (6.23b 
°° 2K, Ce 
Assuming v > 0, the time step for the unsplit scheme for advection-diffusion is, 


2 


(6.23c) 


At ap 


pier 
VAx+2K , 
From which we can see that, 


At,y<At, — (6.24a) 


AD— 
and, 


At,ySAty — (6.24b) 


This leads us to conclude that the corresponding direct (unsplit) scheme will have a smaller allowable time 
step than the split scheme (note also that as Ax gets smaller Atap becomes much smaller than At,). From 
our previous theory of split schemes we can derive a split scheme for the advection diffusion equation 
(5.1). We replace (5.1) by two PDEs, 


U,t+vU, =0  (6.25a) 
and 


U,=K,U,, — (6.25b) 
(6.25a) is the pure advection equation and (6.25b) is the pure diffusion equation. Let La(At,) and Lp(Atp) 
be FD operators for solving (6.25a) and (6.25b) respectively. We will show that Lp(At) La(At) is a (term) 
split operator sequence for solving (5.1) where At = min(At,, Atp). From the standard Taylor expansion 
with x constant, 


2 


Ult + Atx)=U(tx)+At0,U+—0, U+0(At') (6.26) 
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Firstly we use the advection-diffusion equation (5.1) to replace the time derivatives in (6.26) by spatial 
derivatives. (5.1) can be written, 


0,U=(-vo, +K,¢,,JU (6.27a) 
Hence, 
0, U=(0 ,0, )U=(-vo,,+K, 0,, (-v0,,+K, 0,, )U 
=(v’d,,—2vK,@,. +K,’@,,.. U 


x XXX 


(6.27b) 


Substituting (6.27ab) into (6.26) gives, 


JU+O(At?) (6.28) 


X~ XX X~ XXX XXXX 


2 
U(t + At,x)=U(t, x)+At(-vo, + K,0 us W*s,, —2vK,0 +K,"0 


“. U(t + At, x) = Lap(At) U(t, x) + O(At’) ~— (6.29a) 


where Lap(At) is an unsplit differential marching operator to solve the advection-diffusion equation given 
by, 


X 7 XXX xxxx FO(At*) (6.29b) 


At? 
Lap (At)=14+At(—vd, + K,0n.)4 Wx ~2vK 0, +K,’8 
(6.29b) is made into a difference marching operator by replacing all partial derivatives by FD 


approximations. By previous work we know that a differential marching operator, La(At), to solve the 
advection equation (6.25a) is given by, 


2 
L, (At)=I-Atvo , va, +O(At*) (6.30) 


We derive a differential marching operator, Lp(At), to solve the diffusion equation (6.25b). (6.25b) can be 
written, 
0,U=K,0,,U — (6.31a) 


hence, 
0,,U=0,0,U=(K, 0, (K, 0,, JU 


xX XX xX XX 


=K0.U 


XXXX 


(6.31b) 
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Substituting (5.18ab) into (6.26) gives, 


At? 


U(t + At, x}-U(t, x)+AtK,,0,, U+K, 0 U+O(At®) (6.32) 


XXXX 


U(t + At, x) = Lp(At) U(t, x) + O(At’) — (6.33a) 


where Lp(At) is an unsplit differential marching operator to solve the diffusion equation given by, 


2 
Ly (ABH1+AtK,3,.+K,70 +O(At?) — (6.33b) 


XXXX 
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(6.33b) is made into a difference marching operator by replacing all partial derivatives by FD 
approximations. It now remains to show that Lap(At) is approximated by the split sequence Lp(At) La(At). 


2 
Lan (ar) teatk,0, Ki 


XXXX 


2 
104s) [ awe, v0, 10.88) 


2 
=HAtK,0,. Ke 


XXXX x XXX 


2 
Atv, -At?vK 8 v0. +0(At’) 


=L (At) +O(At*) 


Hence the unsplit operator and the split operator sequence are the same up to O(At’). 


Using the previous FD schemes the corresponding split scheme to solve the 1D advection-diffusion 
equation (5.1) is, 


La: u,=u" -Z(ul-ut,) (6.34a) 


—n K Ath ie re 
atl =u; i — 2u,+u:) (6.34b) 


where, 


Output is given in Figure 6.3 for the same conditions as for Figure 5.3. It can be seem that the peak 
concentration is in the correct place as before. This was achieved in fewer time steps. However the shape 
of the concentration profile is slightly different to that given by the unsplit scheme - it is higher and 
narrower. Perfect numerical schemes would give the same results for split and unsplit algorithms but no 
scheme is perfect and the effect of numerical diffusion in the advective solvers plays a different role in 
each algorithm. 
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advection-diffusion: initial profile (dotted) and split numerical solution (+) at later time, 
1 T T, T T T T T T T 


0.9 + 4 


0.8} a 


0.74 t+ | 
0.6 - | 


0.5} 4 


concentration u 


0.4} + | 
0.3} + | 


0.2} 


+ 

+ 

0.1} es 
+ 
= 


0) 40.20.30. 40.«50 S60. 70S BO. 90.-*100 
Figure 6.3 Initial profile and split numerical solution (+) to the 1D advection-diffusion equation at t = 57s. 


6.3.3 Exercise 6b 


1. For v>0 show that, a) lim 
Ax>0 At a Ax0 At ap 


What do you conclude from this? 
2. Go through the derivation of Lap, Lp and L, and show that 
Lap = Lp La + OAC). 


3. Write a program for the split scheme (6.34ab) and reproduce Figure 5.3. 
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7. Systems of Equations 


7.1 Introduction 


The world is governed by natural laws many of which can be expressed as systems of PDEs. An important 
example are the Navier-Stokes equations that, together with the Continuity Equation form a system of 6 
coupled PDEs which describe fluid flow in 3D. These equations are difficult to solve even approximately. 
In the following, a relatively simple system of PDEs has been chosen to illustrate the FD approach. 


7.2 The Shallow Water Equations 


The Shallow Water Equations (SWE) may be expressed in 1D or 2D and provide a simplified model of 
water flow which may be used to simulate many situations including river flow and tsunami propagation. 
For simplicity we consider only 1D. This is still useful and there are many 1D SWE software packages 
used by hydraulic engineers to make flow calculations for real life applications. 


The 1D SWE form a system of two coupled non-linear hyperbolic PDEs with independent variables, 
t (time) and x (distance along the flow) and dependent variables h = h(t,x) (flow depth) and v = v(t,x) 
(flow velocity). 
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Using our usual subscript notation to indicate partial differentiation the 1D SWE can be written 
compactly as, 


EE) 


where, 


F(U)= he. 7. 2b 
—= |B] | hv?+gh?/2? 2?) 


0 
Ur 4s, at 


where g is the acceleration due to gravity, U is the matrix of dependent variables, F(U) is called the flux 
vector and S(U) is called the matrix of source terms which here only consists of bed slope (S,, measured 
positive downwards) and friction (S-) terms. 


7.3 Solving the Shallow Water Equations 


7.3.1 Theoretical Background 
Explicit FD schemes need a limited time step for stability. Previous stability analyses assume linearity but 
the SWE are not linear and so we cannot use these approaches directly. It is therefore necessary to 


examine the SWE in some detail to determine stable time steps before looking at explicit FD schemes. 
Neglecting source terms and using the Chain-Rule (7.1) can be expressed in quasi-linear form as, 


U,tJU,=0 (7.3) 


where J is the Jacobian matrix given by, 


OF, OF, 
oU, 0U 

=| Ee ee. ea) 
6U| | dF, OF, 
aU, dU, 
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It is easily shown that, 
J : (7.5) 
eae +gh 2v 


J has two real eigenvalues, 
A =Vt+ gh ,A,=v-/gh (7.6) 


These eigenvalues are real and distinct and therefore give rise to corresponding linearly independent 
eigenvectors q; and q, (which we write as column vectors). By standard Linear Algebra theory there is a 
transformation of coordinates that diagonalizes J. This is now shown. Let, T = [ qi qp J, then, 


JT=J [ai a] 
=[Jai J q] 


=[Aiqi As gp] (by definition of eigenvalue) 


1.0 


vlasl |), 


0 
= bi 
Or 


therefore, 


W, 
We introduce a change of variable by defining wy | by, 
2 


W=T'U (7.8) 
Assuming that T" is locally constant, then (7.3) becomes, 


TW,+JTW,=0  (7.9a) 
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Multiplying through on the left by T" and using (7.7) gives, 
W,+DW,=0 = (7.9b) 


Since D is a diagonal matrix this system of equations has been decoupled by the change of variable. 
Writing out each row of (7.9b) gives, 


‘=0 = (7.10a) 


= +(v+/gh ) 


OW, 
0 0 


Xx 


0 


W, OW. 
+(v—,/gh)—-=0 (7.10b 
es (v-/g as ( ) 


(7.10a) and (7.10b) are two linear advection equations with two wave speeds (v+,/gh),(v—,/gh ) 


respectively. Any explicit scheme to solve the SWE must take account of the above wave speeds to ensure 
stability. Note that these wave speeds vary over both space and time. 
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7.3.2 Heuristic Time Step Calculation 
The following is a treatment of the time step calculation for the SWE based on the previous heuristic 
analysis of Appendix C. Suppose that a given explicit scheme to solve the advection equation u, + v ux = 0 


is stable when, 


At<Ax/(5v) (7.11) 


where v is the wave speed (assume v > 0) and At and Ax are the time and space steps respectively. By the 
previous analysis, the same scheme will be stable when solving the SWE if the maximum wave speed 
(chosen over the spatial interval at a given time level) is used in the stability inequality. 


ie. At< 


(7.12) 


Ax 
Smax(v+ Jgh)(v—eh)} 


Furthermore, because the wave speeds also vary with time, this constraint must be applied at each 
time step. 


7.4 Example Scheme to Solve the SWE 


The SWE (7.1) with S(U)=0 look very similar to the linear advection equation, u,; + v u, = 0. We see that 
U, is replaced by u, and F, is replaced by v u,. The FOU scheme (Chapter 4) to solve the linear advection 
equation is, 


n n 
n+l __ n_Atv(u; —u;,) 
j TU; a 


(7.13a) 


The FOU scheme can be shown to be stable (for v>0) if, 


At<Ax/v.  (7.13b) 


We apply the same scheme (in the sense that the time derivative is replaced by the first order forward 
difference approximation and the spatial derivative is replaced by the first order backward difference 
approximation) to solve the SWE giving, 


Xx 


1 —l 
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By the previous analysis the time step inequality (7.13b) becomes, 


A (7.14b) 


Ax 
t< 
max{(v-+ Jgh)(v—Jgh)} 


(7.14b) must be applied at each time step. 


The above analysis shows that a program written to solve the linear advection equation may be modified 
to solve the SWE by appending an extra row to each array variable matrix. Prior to each time step a 
routine to calculate the new time step must be invoked. After each time step h and v are found from U 
(h=U;, v= U2/U;) and then U and F are initialized before the next time step. 


Figure 7.1 shows the Lax-Friedrichs scheme solution to the SWE for a collapsing water column (not 
shown is the graph for the associated water velocity). Initially the water is still and the surface profile is 
given in Figure 7.1. The program was run for 6 seconds using 201 grid points over [0, 200] and a time step 
safety factor of 0.95. Zero-gradient boundary conditions were used for both water height and velocity. 
Note that this problem has discontinuous initial conditions for h that can cause classical schemes problems 
(see the oscillations in Figure 7.1). There exist modern schemes to cope with discontinuities but they are 
beyond the scope of this text. 


127 


Figure 7.1 Lax-Friedrichs solution to the Shallow Water Equations for a collapsing water column. Initial 
conditions (__), numerical solution (---) at t = 6s. 
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Code to implement the above scheme may be downloaded. A Case Study to illustrate the use of the 1D 
SWE is given on the website. This study also gives useful information on how to set up a problem and 
assess results. 


7.5 Exercise 7 
1. Show that J is as given and use the Chain Rule to express (7.1) in the form of (7.3). 
2. Find the eigenvalues of J. 


3. Given that an explicit scheme to solve the linear advection equation with wave speed v is stable for 
At<2Ax /(3v) state the (heuristic) stability criteria when this scheme is used to solve the SWE. 


4. Comment on the sign of the wave speeds for the SWE when the flow is such that v > (gh)'”. 


5. Modify one of your existing codes to solve the linear advection equation using the Lax-Friedrichs 
algorithm (use previous initial and boundary conditions). Verify your code. 


6. Modify your program in Q5 to solve the SWE using the conditions of Figure 7.1. Reproduce Figure 
dls 


7. Animate the solutions for height and velocity in Q6 for and run until the water passes out of the 
domain. 


8. Repeat Q7 with solid left and right hand boundary conditions (see Appendix B). 
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Appendix A: Definition and Properties of Order 


A.1 Definition of O(h") 


For our purposes, 
f(h) = O(h") 


(pronounced, ‘f of h is order h to the n’) means, 


limt)_ C, (Al) 


h>0 h 


where C is a non-zero constant. 
500h°+3h*-2h _ 
h 


e.g. 500 h° + 3h* - 2 h =O(h) because, lim —2 
> 


4 


9h 
4 = 4 i _————— 
e.g. 9h” = Oh") because, lim nf 9 
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A.2 The Meaning of O(h") 


If f(h) = O(h") then, for small h, equation (A.1) gives, 


EO) ag. 
h" 


f(h) = Ch® (A.2) 
Equation (A.2) says that for small h, an error which is O(h") is proportional to h”. In particular if the error 


is O(h) then it is proportional to h which means that halving h halves the error. If the error is O(h’) then it 
is proportional to h? which means that halving h reduces the error by a factor of 27 = 4. 


A.3 Properties of O(h") 


In the analysis of finite difference schemes we will need to use the following properties of Order notation. 
Let f(h) = O(h"), g(h) = O(h™) where 0 < m < n. Let K be a non-zero constant. Then, 


A.3.1. K f(h) = O(h") 
A.3.2. f(h) + g(h) = O(h”) 
A.3.3. f(h)/h = O(h™') 
A.3.4. f(h)/g(h) = O(h""™) 
A.3.5. f(t) g(h) = O(h"™) 
Proof of A.3.3 


By definition of Order we have to show that, 


. f(h) : 
li =— is a non-zero constant. 
h>0 =f" 
Since f(h) = O(h"), 
. f(b 
lim m = C, where C is a non-zero constant. 
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Se pe 


“450 he hoo ” 


=C 


End of proof of A.3.3 


A.4 Explanation of the Properties of O(h") 


A.3.1 says that scaling a function by a constant doesn’t change its order. In particular f(h) and —f(h) have 
the same order. 


A.3.2 says that the order of the sum of two functions of different orders is the smaller of the orders of the 
two functions. e.g. 


O(h’) + O(h*) = O(h’). 


A.3.3 says that dividing a function by h reduces its order by 1. A.3.3 is a special case of A.3.4 that says 
that dividing a function by a function of order m reduces its order by m. e.g. O(h®)/ h? = O(h’). 


A.3.5 says that the order of the product of two functions is the sum of their orders. e.g. 
O(h*) O(h?) = O(h?). 

A.5 Exercise A 

1. Find the order of the following functions: 

2. a)2x+5x*>-3x° b)4x’-17x*+3x* oc) sin(x) d)x 


3. Prove A.3.1) — A.3.5) 


Download free ebooks at bookboon.com 


114 


Introductory Finite Difference Methods for PDEs Appendix B: Boundary Conditions 


Appendix B: Boundary Conditions 


B.1 Introduction 


When solving a PDE using a finite difference (FD) scheme we may need to specify ghost grid points and 
associated ghost values for the dependent variable at these points. 


e.g. in the FOU scheme to solve the 1D linear advection equation (Chapter 4) we need a ghost point to the 
left of the first grid point. This is illustrated in the following diagram: 


ghost point 
Ax 
<<» 
Ax 
<0} i 
i=N 
i=0 i=l i=2 i=3 


Figure B.1 Ghost (white) and grid (black) points for the FOU scheme 


Notes: 


1. In general for a 1D region on which grid points are indexed by i = 1, 2, ... , N, we will index a left 
ghost point by i = 0 and a right ghost point byi=N +1. 


2. In many computer languages (e.g. Scilab and Matlab) array indices start at 1 so we must shift indices 
when using them to code up a FD scheme with a left ghost point. For example to code the FOU 
scheme we shift grid indices so that the left hand ghost point index is | and so indices 2 to N+1 
represent the N computational grid points. If a scheme (e.g. Lax-Friedrichs) requires both left and 
right ghost points then the left ghost point has index 1 and the right ghost point has index N+2. 
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B.2 Boundary Conditions 


In a computational region (which may be 1, 2 or 3D) ghost points and associated values occur at or 
adjacent to the boundaries of the region. Conditions leading to the prescription of ghost values are called 
boundary conditions. Boundary conditions are derived from the underlying physics of the situation. The 
correct treatment of boundary conditions is vital for accurate problem simulation. This is illustrated by 
considering the grid for a 2D region shown in Figure B.2. Grid points are indexed by i= 1, 2, ... , N in the 
x direction and by j = 1, 2, ... , M in the y direction. Indices 0, N+1 and M+1 indicate ghost points. To see 
the effect of different boundary conditions we consider the region to represent two different problems 
(modelled by the same PDEs who solution is approximated by a FD scheme requiring the ghost points 
shown). 
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B.2.1 Scenario 1 


Consider the region in Figure B.2 to represent a channel of water viewed from above and flowing from 
left to right. The rows having indices j = 1 and j] = M correspond to the right and left sides of the channel 
respectively. If water cannot flow through or over the sides of the channel we must impose solid boundary 
conditions at each side by specifying the values of water depth and velocity at these locations in a special 
way. The columns with indices i= 1 and i=N are at inflow and outflow boundaries respectively and 
require water depths and velocities to be specified according to flow type. 


B.2.2 Scenario 2 


Consider the region in Figure B.2 to be a near-shore area of ocean as viewed from above with waves 
travelling from left to right and impacting on a solid harbour wall that cannot be overtopped. The harbour 
wall, being solid, requires a solid boundary condition be imposed on the column with i = N. The column 
with 1 = 1 requires a time dependent boundary condition (for water depth and velocity) to generate the 
incoming waves. The rows with j = | andj = M are edges of the finite computational domain and require 
transmissive boundary conditions to be imposed there so that waves can pass in and out of the 
computational domain. 


(0, M+1) (N+1, M+1) 
- ots oo = 


(1, 0) (2, 0) (3, 0) (N,0)  (N+1, 0) 


Figure B.2 Ghost (white) and grid (black) points for a 2D region 
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From the two previous scenarios it is clear that changing the boundary conditions changes the problem. 


B.3 Specifying Ghost and Boundary Values 
There are two basic ways to specify values of the dependent variable at ghost points. For clarity we will 
assume that the ghost point is at index i= 0 in a 1D domain and hence i = 1 is the index of the grid point 


on left hand boundary. 


B.3.1 Dirichlet Boundary Conditions 


In Dirichlet boundary conditions the value, ug , of the dependent variable at a ghost point is specified in 


some way. Examples include: 


a) U,j=constant.e.g. uj =0 in the FOU scheme for the 1D linear advection equation (this condition 


indicates that there is no more pollutant entering the river, see Chapter 4). 


b) uj=f(n), which is a time dependent boundary condition (e.g. a tidal boundary as in the previous 


harbour example). 


c) Uj j=U,,. This is a periodic boundary condition. This means that what passes out of the right 


boundary will pass in from the left boundary as though the two boundaries were joined together. 
B.3.2 Derivative (von Neumann) Boundary Conditions 


As the name suggests derivative boundary conditions specify the rate of change of the dependent variable 
at the grid point adjacent to the ghost point (i.e. at i= 1 in our case). This can be done in two ways: 


B.3.2.1 U,=f(U). 


Here the derivative of U in the x direction is specified at the boundary grid point i= 1. From this 


information ug can be calculated. One way is to estimate U, at i= 1 by a central difference giving, 


n n 
U,—Up 


U,=f(U)= 
x=f(U) a 


(B. 1a) 
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which rearranges to give, 


uj=u,—2Axf(U) (B.1b) 


B.3.2.2 U,=f(U). 


Here the derivative of U in the direction of n, the outward pointing normal to the boundary is given at the 
grid point adjacent to the ghost point. Note that this direction is opposite to the x direction at the left hand 


boundary (i= 1). From this information uj can be calculated. As before we will use the left hand 


boundary and estimate U,, at i= 1 by a central difference giving, 


n n 
Up—U, 


U_xf(U)=—+ 
en) 2Ax 


(B.2a) 


which rearranges to give, 


uj=u5+2Axf(U) (B.2b) 


Notes: 


1. (B.1b) and (B.2b) are different because the directions of the given derivatives are opposite (the 
formulae would be the same at the right hand boundary). It is more usual to use (B.2b) then (B. 1b). 


2. Inthe derivative boundary condition examples central differences were used but other estimates 
could easily have been used (e.g. first order backward difference for the left hand ghost value). 


3. The formal accuracy of a scheme may be reduced if ghost values are calculated on the basis of a 
method whose accuracy is less than the spatial accuracy of the scheme. e.g. a spatially third order 
scheme may drop to second order if a central difference (second order) is used to calculate a ghost 
value. 


4. In many numerical simulations both Dirichlet and von Neumann boundary conditions are used. These 
are called mixed boundary conditions. 
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B.4 Common Boundary Conditions 
We give some boundary conditions that are used frequently in problems. 
B.4.1 Transmissive Boundary Conditions 


Computational domains are finite and it may be that we need quantities simply to pass out of a boundary 
when they reach one that is not solid. This is often done by specifying a zero gradient normal to the 
boundary in the variable of interest on the boundary i.e. we use a derivative boundary condition in which 
the derivative is zero. As an example consider Figure B.2 with Scenario 2 in which water height h is a 
dependent variable and the solver requires a ghost value at i = 3, j = 0. Grid point i=3,j=1isona 
transmissive boundary so we want waves to pass through. To impose the transmissive boundary condition 
on water height, h, we set h,=0 at grid point i = 3, j = 1. Using a first order backward difference 
approximation, 

oy h3 hs, 


h,=0* B.3 
: ey (B.3a) 


“ht, =h2, — (B.3b) 


What do you see? 


A skateboarder? 


We see an All Options Junior Trader leaving 
the office after a successful day. 
He has the freedom to make his own decisions. 
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B.4.2 Solid (reflective) Boundary Conditions 


There is no flow through a solid boundary. This implies that the normal component of flow velocity at a 
solid boundary is zero. This condition is implemented by copying the normal component of velocity at the 
neighbouring interior grid point to the ghost cell and changing its sign. As an example consider Figure B.2 
with Scenario 1. Grid point i=3, j=1 is on a solid boundary and we wish to find the velocity at the ghost 
point i=3, j=0. Suppose the velocity at neighbouring interior grid point i=3, j=2 has x and y components vx 
and v, respectively. The component of velocity normal to the solid boundary is v, so the velocity at the 
ghost point i=3, j=0 has x and y components v, and -vy respectively. Another way of obtaining this result 
is to linearly extrapolate using the normal velocities at i=3, j=2 (i.e. vy) and i=3, j=1 (i.e. 0). Variables 
other than velocity at ghost points next to solid boundaries are usually found by a zero gradient approach. 


B.4.3 Slip and No-Slip Boundary Conditions 
At a solid boundary we have a choice of tangential velocity component. If there is appreciable friction at 
the boundary then it is said to be a no-slip boundary and the tangential velocity component is zero (along 


with the normal component). If friction is not present then the boundary is said to be a slip boundary and 
the tangential component of velocity may be extrapolated from interior tangential components. 


B.5 Exercise B 


1. A 1D PDE is to be solved numerically. The computational domain consists of N grid points 
numbered 1, 2, ... , N. In each case give the indices of the ghost grid points for the following FD 
schemes (which can be found in Chapter 4), 


a. FTCS scheme, b) Crank-Nicolson, c) Lax-Wendroff, d) Lax-Friedrichs. 

2. Run the FOU scheme for the 1D linear advection equation (Chapter 4) using a periodic boundary 
condition and describe how the solution evolves in time. (Remember indices start at 1 in Scilab or 
Matlab). 

3. The computational domain [0, 100] is discretised using 101 points. Initial values are given by 


U(0, x) = sin(x). The FOU scheme is used. Calculate the ghost point and the initial associated u value 
given the following boundary conditions, 


a) Periodic, b) Dirichlet: uj=0.1Lu,,,,=—0.1 c) Derivative: U,=0, 


d) Derivative: U,,=0.1. e) Derivative: U,=U. 
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b) 


c) 


Repeat Q3 for the Lax-Friedrichs scheme. 

In Scenario | the flow velocity vector at grid point i=2, j=M-1 is (vx, Vy) = (4, 8). 
Calculate the flow velocity vector at ghost point i=2, j=M+1. 

Assuming a no-slip boundary find the flow velocity vector at grid point 1=2, j=M. 
Repeat b) assuming a slip boundary. 


In Scenario 2 the water height at grid point i=3, j=1 is 8. Find the water height at ghost point i=3, j=0. 
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Appendix C: Consistency, Convergence and Stability 


C.1 Introduction 


Consistency, convergence and stability are important concepts but much of the theory has only been 
developed for special cases so a combination of theory and numerical experimentation is often the only 
way to proceed (e.g. to find the maximum stable time step). For clarity we restrict our attention to a single 
PDE where there are two independent variables, t and x. The dependent variable is U = U(t, x). In operator 
notation the PDE is written, 


LU=0 (C1) 


Lis a differential operator and U is the analytical (i.e. exact) solution of (C.1). 


In order to approximate the solution of (C.1) the computational region is discretised into a finite set of grid 
points, x;, where Ax is the (constant) grid spacing. The approximate solution is found at a set of time levels, 
t" where At is the (variable) spacing between time levels. The analytical (i.e. exact) solution to (C.1) at 

(t", x;) namely U(t", x;), is denoted by U;”. In operator notation the Finite Difference (FD) scheme to 
approximate the solution of (C.1) is written, 


Darax¥; =9  (C.2) 


D is a difference operator and u;" is the solution to the FD scheme at (t", xj). 
Notes: 


1. The idea of the FD scheme (C.2) is that u;" approximates U;' and the approximation becomes better 
and better as Ax and At become smaller. Let, 


ei = u; - U; (C.3) 
e;' is called the pointwise error (also called the ‘global error’). 


2. Initially (i.e. time level 0) U is known at all grid points and u;,’is taken to be U;’ so e;° = 0 at all grid 
points. As iterations of the FD scheme introduce errors, in general e; #0. It may be that as iterations 


continue errors are compounded and e;' grows unboundedly making the FD scheme useless. 
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The following concepts look at properties of the FD scheme with respect to errors. 


C.2 Convergence 


The FD scheme (C.2) for the PDE (C.1) is said to be convergent at (t, x) if the pointwise error at (t, x) 
tends to zero as Ax and At tend to zero, i.e. 


e; >0 as (t",x,) (t,x) when At,Ax—>0. (C4) 


To be useful our FD scheme must be convergent but this is very difficult to prove for many schemes. 
Convergence implies that a solution of the FD scheme approximates a solution of the PDE. 


C.3 Consistency and Scheme Order 


A measure of how well the exact solution of the PDE satisfies the FD scheme is given by the 
truncation error, 


Digg? CS) 


The FD scheme (C.2) is consistent with the PDE (C.1) if the truncation error tends to zero as Ax and At 
tend to zero, 1.e. 


DaraxU? 0 as At, Ax—>0 at (t",x;). (C.6) 


Consistency is necessary for convergence. 

The order of the truncation error is obtained by Taylor expansion of its terms about (t", x;). In many cases 
the order of the truncation error is the same as the order of the pointwise error so that the truncation error 
is a good (and accessible) guide to the accuracy of a FD scheme. We define the formal order of accuracy 
of a FD scheme by the order of its truncation error. 


C.3.1 Example Calculation of Consistency and Scheme Order 


The 1D linear advection equation is, 


U,+vU,=0 (C7) 
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The FTCS scheme to solve (C.7) is given in Chapter 4. Writing this scheme in the form of (C.2) gives, 


a; =U, +V Ui —U; = 0 (C.8) 
At 2 Ax 


To determine scheme consistency and order we replace u in (C.8) by the exact solution, U, of (C.7) to give 
the truncation error, 


ye =0 (C9) 
At 2 Ax 


Terms in (C.9) are replaced by their Taylor expansion about the i” spatial point and the n" time level to give, 


U? +At U,+O(At?) -U? i 


2 
(| u'+4axu, 48 U, +0(Ax?) 
At 2Ax 2 


(C.10) 
n Ax?* 3 
= D. es ) ) 


Simplifying (C.10) (using properties of O notation) gives, 


U, + v U, + O(Dt) + O(Dx’)  (C.11) 
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and using (C.7) the truncation error is finally, 


O(At) +O(Ax’) — (C.12) 


Clearly the truncation error (C.12) tends to zero as At, Ax +0 which demonstrates that the FTCS scheme 
(C.8) is consistent with the original PDE (C.7). Furthermore the truncation error is first order in time and 
second order in space so the FTCS scheme is said to be (formally) first order in time and second order in 
space. 


Notes: 


1. Inthe Taylor expansions we expanded to second order for t and to third order for x. It doesn’t matter 
if we expand too far as additive order expressions collapse to the term of lowest order. 


2. Consistency does not imply that a FD scheme is any good! In fact the FTCS doesn’t work for reasons 
we now explore. 


C.4 Stability 


We have seen from numerical results in Chapter 4 that care must be used when choosing the time step, At. 
A scheme that produces acceptable results for small At can give results that grow unboundedly as 
iterations continue if too large a At is chosen. In this case the scheme is said to have become unstable. A 
FD scheme is stable if and only if pointwise errors do not grow unboundedly with time which will be the 
case if, after some time level, N, 


n 


e 


,n2=N. (C.13a) 


Another way of defining stability comes from the following argument. Rewriting (C.3) gives, 


ur = U; a e (C. 13b) 


For stability e;" must not grow unboundedly and since U;’ is finite, (C.13b) implies that u;" must also not 
grow unboundedly which is true if, 


n 


u 


»n=N. (C.13c) 


The Lax Equivalence Theorem says that for a /inear PDE, a consistent FD scheme to approximate its 
solution is convergent if and only if it is stable. 
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C.4.1 Heuristic Stability Analysis 


The following reasoning gives a guide to choosing a stable time step and may be applied to PDEs where 
there is a speed of propagation. We apply our reasoning to the simple 1D linear advection equation of 
Chapter 4 (we will see in Chapter 7 that this reasoning extends to more complicated PDEs). Reasoning 
goes as follows: during time step, At, the concentration profile travels a distance |v|At downstream. To 
capture the solution numerically it seems reasonable not to let the concentration profile move more than 


one grid interval, Ax, in a single time step, At. This means that, | v |At<Ax , 


ame. (614g 
|v | 


Inequality (C.14a) gives a heuristic guide to the maximum allowable time step (for the linear advection 
equation). It is often a good idea to multiply this maximum value by a ‘safety factor’ F where, F < 1. 
Numerical experiments can then be used to determine F. It is customary to let, 


At 
c= (C:14b) 
Ax 


c is called the Courant number. As we will see, schemes are often stable for c less than or equal to some 
number. (C.14a) implies that, 


|c|Sl  (C.14c) 
This is called the CFL condition (after Courant, Friedrichs and Lewy). 
C.4.2 Matrix Stability Analysis 
This is a more mathematical approach to stability analysis but it only applies to /inear FD schemes. By 


writing out a linear FD scheme at each grid point and expressing the resulting set of linear equations as a 
single matrix equation, a linear FD scheme with 2 times levels can be written as, 


Au™'=Bu" (C15) 
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where, u"=(u? A silly y and A and B are NxN matrices (A =I for an explicit scheme). If the FD 


scheme is consistent it is satisfied by the exact solution of the PDE (neglecting the vanishing truncation 
error) and so, 


AU;"=BU; (C16) 


Subtracting (C.16) from (C.15) and using (C.3) gives, 


Ae™'=Be" (C.17a) 


se = A'Be" (C.17b) 


| ntl 


=|A "Be" 


(C.17¢) 


| ntl n 
“Je e 


(C.17d) 
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Hence by (C.13a) the FD scheme (C.15) is stable if, 


A 'B| <1 (C18) 


The matrix norm in (C.18) is induced by the vector norm in (C.13a). Hence the stability of a (linear) FD 
scheme can be investigated by finding the norm of the matrix AB. This is the matrix method for stability 
and may be quite difficult to implement. It should be noted that there are many definitions of norms and a 
FD scheme may be stable in one norm but not in another. 


C.4.2.1 Example of Matrix Stability Analysis 


From Chapter 4 the FOU scheme to solve the 1D linear advection equation (C.7) can be written, 


“t=cu",+(1-c)u? (C.19) 


1 


u 


Where c is defined in (C.14b). We use matrix stability analysis to find the maximum time step for scheme 
stability. Writing (C.19) out for each grid point, assuming u; = 0 at all time levels and expressing the 


resulting linear equations in the form of (C.15) gives, 


(l-c) 0 0 .. 0 
c (l-c) 0 

u's) 0 c (l-c) 0 = 0 |u" (C.20) 
0 c  (I-c) 


(C.19) is stable if the norm of the above matrix is less than or equal to 1. If we use the infinity norm for 
vectors (which is the maximum absolute value of components) the induced matrix norm is the maximum 
of the sum of all absolute values in each row. For our matrix this is |c| + |1-c|. For 0 < c<1 (which implies 


that v>0), |c| + |1-c| =c + (1-c) = 1. Hence the FOU scheme is stable for 0 < c<1 and hence the maximum 
allowable time step for stability is = which agrees with the previous heuristically derived result. The 

Vv 
FOU scheme is said to be conditionally stable. 


C.4.3 Von Neumann Stability Analysis 


This analysis of stability is due to von Neumann and is based on (C.13c). We assume that the FD scheme 
is linear and that boundary conditions are periodic. Using the complex version of Fourier’s Theorem and 
rescaling so that u has period 27, we may write, 
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i= Dae (C.21) 


k=—00 


where, 


j=v-—l, and ler is the amplitude of the k" Fourier component. 


By linearity it is enough to examine the behaviour of a single Fourier component so we replace 


u’ bygte™™ in (C.13c) and rearrange to get, 


n+l 


Bx 


k 


G=2*<1_  (€.22) 


This is the von Neumann condition for stability. G is called the amplification factor. Von Neumann 
stability is carried out by the following steps: 


Step 1. Replace each instance of U,’ in the FD scheme by its corresponding single Fourier component. 
Step 2. Rearrange to get G. 


Step 3. Use the constraint (C.22) to obtain the condition for At (this step could be algebraically tricky). If 
(C.22) can never be satisfied for At > 0 the scheme is unconditionally unstable. 


C.4.3.1 Von Neumann Stability Analysis: Example | 


The FTCS scheme for the 1D linear advection equation is, 


n+ n vAt n n 
uo ek 2 ue), (C.23a) 


which can be written compactly as, 


n+l 


u=u}—S(u-ul), C230) 


where c is the usual Courant number defined by (C.14b). 


Step 1: The FD scheme is linear and we assume periodic boundary conditions. Replacing each term in 
(C.23b) by its k" Fourier component gives, 
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n+l jkx; nikx; © (on kx: ny jkx\_) 
gitleii tei —S{ete™ —grel ) (C.24) 
Step 2: Noting that x;.) = x; + Ax, x;.) = x; - Ax, gives, 


ntl. jkx; nn . jkx; Cin jk(x;+Ax) n_, jk(x;-Ax) 
& © =B,e -S{gte —8,& ) (C.25) 


Dividing through by gre" gives, 


n+l 

BS. 4 Of jkax —jkax 

2s _=]-—|e"™™ -e (C.26) 
gy 3 
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Step 3: From (C.22) for stability we must have, 


Ge= eee }a (C.27) 


jO_ ,-j0 
Using the well known identity, sino)4{ S™ 
J 
(C.27) gives, 


|I-jCsin(kAx)|s1_ (C.28) 


Squaring each side and evaluating the squared modulus gives, 


1l+c’sin?(kAx)<1  (C.29) 


This inequality can only be satisfied (for all k) if c = 0 which implies that At = 0. i.e. there is no feasible 
value for the time step that makes the scheme stable. i.e. the FTCS scheme for the 1D linear advection 
equation is unconditionally unstable. The FTCS scheme is therefore useless even though we have shown 
that it is consistent! 


C.4.3.2 Von Neumann Stability Analysis: Example 2 
We apply von Neumann stability analysis to the FOU scheme (C.19). 


Step 1: Replacing each term in (C.19) by its k™ Fourier component gives, 


pte =(1-c)gte™ +cgrei (C.30) 


Noting that x;_; = x; - Ax, gives, 


grttelis —(J-c)gte™™ +egre™™) == (C31) 
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Step 2: Dividing through by g'e"™ gives, 


Sk -(]-c)+ce™* (C.32) 
Step 3: For stability we must have, 


G = |(-c)+ce™ 


<1 (C.33) 


The well-known triangle inequality states that, |a+b |S| a |+| b], 


hence, 


\(l-c)+ce™ 


<|(1-c)}+}ce*™ 


=((1-c)|+|c] (C34) 


When O<c<I,| 1—c|=l—c and | c |=c, therefore (C.34) gives, 


\-c) +ceu™ 


<|(1-c)}+|c|=(I-c)+e=1 — (C.35) 


Hence the FOU scheme for the 1D linear advection equation is stable when 0<c<1 which means that 


A 
At< = . The FOU scheme is said to be conditionally stable. 
Vv 


Notes: 
1. Stability analysis hasn’t been worked out for most non-linear schemes (PDEs). 


2. Strictly speaking, the von Neumann stability analysis requires periodic boundary conditions but 
seems to work even when this is not the case. 


C.5 Exercise C 


1. Show that the FOU scheme for the 1D linear advection equation is consistent and find its formal 
order. 


2. Repeat Q1 for a) Lax-Friedrichs b) Lax-Wendroff c) Crank-Nicolson 
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3. Show that the 5-point scheme of Chapter 3 for the 2D Laplace’s Equation is consistent and find its 
formal order. 


4. Using heuristic analysis, estimate the maximum allowable time step for an explicit scheme to solve 
the 1D linear advection equation for pollution in a river 0.5 km long flowing at 4m/s using 200 grid 
points. 

5. Use von Neumann stability analysis to show that the FD scheme to solve the linear advection 
equation for v > 0 using first order forward differences in both space and time is unconditionally 


unstable. 


6. Use von Neumann stability analysis to show that the Crank-Nicolson scheme to solve the 1D linear 
advection equation is unconditionally stable. 


7. Use von Neumann stability analysis to investigate the stability of the Lax-Friedrichs scheme to solve 
the 1D linear advection equation. 


8. Repeat Q7 for the Leap-Frog scheme. 


9. Repeat Q5-Q7 using matrix stability analysis. 


6007 'A'G Swewsds Wy] 222] @ 


iB \ 
3 It’s only an 


opportunity if 
you act on it 


IKEA.SE/STUDENT 


eel es aaa 


Download free ebooks at bookboon.com 


134 


Introductory Finite Difference Methods for PDEs Appendix D: Convergence Analysis for Iterative Methods 


Appendix D: Convergence Analysis for Iterative 
Methods 


D.1 Introduction 


Iterative schemes for matrix inversion do not necessarily converge so we need to determine conditions for 
convergence. We would also like to know how fast they converge. In the following we assume that the 
reader is familiar with some standard results from Linear Algebra. We are trying to solve the system of 
linear equations (refer back to Equations (3.6), (3.9)), 


Au=b (D1) 


where A is an NxN matrix, u is a column vector of N unknowns and b is a column vector of N known 
constants. It is always possible to scale each equation so that every entry on the main diagonal of A is 1. A 
can then be written as, 


A=-L+I-U (D2) 


where I is the NxN identity matrix and L and U are NxN lower and upper triangular matrices respectively. 
Substituting (D.2) into (D.1) and rearranging gives, 


u=(L+U)u+b (D3) 
Equation (D.3) is the basis for the following analysis. 
Notes: 
1. Research into computationally efficient ways to invert a matrix continues. Computational efficiency 
is not simply a matter of reducing the number of calculations for a particular class of problem; it also 
depends on computer architecture. It may be that some less sophisticated method is faster than a more 


modern method when run on a parallel computer. 


2. There are many efficient freely downloadable matrix inversion programs so it is almost never worth 
writing your own. 


3. There is probably no best solution to computational efficiency. Where speed is an issue it will pay to 
experiment with several methods and tune them (if necessary) by numerical experiments on small 
problems. 
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D.2 Jacobi Iteration 
Equation (D.3) suggests the iterative scheme, 


ao = (L aie U) rT Ae b (D.4) 


which is the Jacobi iteration scheme in matrix form. (L + U) is called the iteration matrix. In order to 
analyse convergence of the Jacobi iteration scheme we need to look at how the error behaves between 
iterations. Clearly we want the error to decrease to zero as iterations continue. The exact solution to (D.1) 
is u. Let the error after the m" iteration be e”, SO, 


e"=u"—u_ (D5) 


Note that e™ is a vector in RN whose components are the errors at each grid point after the m" iteration. 
Subtracting (D.3) from (D.4) gives, 


u™!-u=(L+U)u"+b-(L+U)u-b (D.6) 


ce™l=(L+U)e™ (D7) 


writing e° for the initial error in the initial (guessed) values for u’, after 1 iteration (D.7) gives, 


e'=(L+U)e’. 
A second iteration gives, 
¢=(L+U)e! 
=(L+Uye’. 


So after n iterations we have 


e"=(L+U)e’ (D8) 


For convergence of the iterative scheme all components of e" must approach zero in the limit, 


i.e. lime” =0 


n> 


~. lim(L+U)"e° =0 


no 
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Since e° is a non-zero constant vector we must focus our attention on the behaviour of (L + U)". We 
assume that (L + U) has n linearly independent eigenvectors, vj), Vo, ... , Vw with corresponding 
eigenvalues, A), Ao, ... , An. By a standard result from Linear Algebra (L + U)” has the same eigenvectors 
with corresponding eigenvalues, 4", Ax", ... , An". The set of eigenvectors form a basis for R‘ so for some 
constants a; we may write, 


0 
e€ =a, Vi +a V2 +... + an Yn, 


. (L+ Ue’ =(L+ UY (ai vitave +... + an Yn) 


=(L+U)"ajvit+(L+U)yawt...+(L+U)" ay w 


=a, (L+U)'vyjta(L+U) wt... + ay (L+U)' w 


=a) Ay" Vi + a2 Ag" Wo +... + an An’ VN 


which will clearly tend to the zero vector if and only if, 


Ai <1, fori=1,2,...,N. 


Definition: The dominant eigenvalue of a matrix is the eigenvalue with the largest modulus. 

Hence we can say that the Jacobi iterative scheme converges if and only if the dominant eigenvalue of its 
iteration matrix has modulus less than 1. We denote the dominant eigenvalue of the Jacobi iteration matrix 
by 2. 


D.3 Gauss-Seidel Iteration 


Using the previous notation it can be shown that the Gauss-Seidel iterative method can be expressed as, 


y™! = U uw L y™! b (D.9) 
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and a similar analysis to the Jacobi iterative scheme shows that, 


er a U e™ +L em (D.10a) 


ve™l=(1-Ly'Ue™ (D.10b) 


ee™lS(I-Ly'Uj™e® (D.10c) 


(I—L)' U is called the Gauss-Seidel iteration matrix. By an exactly similar analysis to that for Jacobi 
iteration, the Gauss-Seidel scheme converges if and only if the dominant eigenvalue of its iteration matrix 
has modulus less than 1. It can be shown that the dominant eigenvalue =. 
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D.4 SoR Iterative Scheme 
A similar analysis to the above shows that the SoR iteration matrix is, 


(I-wU)'(1-—w)I+wL) (D.11a) 


It can be shown that its dominant eigenvalue, is, 


2 


144/1-2 


As before the scheme converges if and only if this has modulus less than 1. 


-1  (D.11b) 


D.4.1 A Special Case for SoR 
For a rectangular pAx by qAy computational region the optimal value of the SoR relaxation parameter can 
be shown to be, 


w (D.12a) 


2 
ae ee ae 
where the dominant eigenvalue, p, of the corresponding Gauss-Seidel scheme is, 
2 
cos scos) 
p 


q 
= D.12b 
My 4 ( ) 


D.5 Theory for Dominant Eigenvalues 


Convergence of iterative schemes depends on the dominant eigenvalue of the associated iteration matrix 
having a modulus less than 1. In general it is difficult and/or computationally expensive to find 
eigenvalues. The following theorem is a quick way to find an upper bound for the modulus of the 
dominant eigenvalue of a matrix. 


D.5.1 Gershgorin’s Theorem 


The modulus of the dominant eigenvalue of a matrix is less than or equal to the sum of the modulii of the 
entries in any row or column. 


03 -01 O1 
eg.Let A=|}0.2 02 0.3 
04 0.1 -0.3 


Download free ebooks at bookboon.com 


139 


Introductory Finite Difference Methods for PDEs Appendix D: Convergence Analysis for Iterative Methods 


The sums of the modulii are: 

Row 1: |0.3|+|-0.1]+/0.1| = 0.5, Row 2: |0.2|+|0.2|+|0.3| = 0.7 
Row 3: |0.4|+|0.1| |-0.3] = 0.8, Col 1: |0.3|+|0.2|+|0.4| = 0.9 
Col 2: |-0.1|+|0.2|+]0.1| = 0.4, Col 3: |0.1|+|0.3|+|-0.3| = 0.7 


Hence the dominant eigenvalue of A is less than or equal to 0.9. If A were an iteration matrix then the 
iteration scheme would converge. If the maximum value of the sum of the modulii of the row or column 
elements is greater than | then the theorem is no use for determining whether the scheme converges. In 
this case we need a way of estimating the dominant eigenvalue. 


D.5.2 Power Method for Estimating Dominant Eigenvalues 


This is an efficient way of estimating the dominant eigenvalue (and associated eigenvector) of a matrix A. 
Start with an arbitrary non-zero vector v’ and define the iterative scheme, 


vt=Avi (D.13) 


It can be shown that as the iteration index i tends to infinity, v'"! tends to A v' where A is the dominant 
eigenvalue of A with associated eigenvector v' and where after each iteration the resulting vector, vil is 
scaled by a constant k; (which is the reciprocal of its first entry) so that its first component becomes 1. The 
i+] 


distance between the scaled versions of v'"' and v' is found and the iteration stops when this distance is 


less than some predefined tolerance. The resulting estimate for the dominant eigenvalue of A is k;. 


D.5.2.1 Example Power Method Calculation 


21 
e.g. af ; . We use the power method with a tolerance of tol = 0.01 to find the dominant eigenvalue 


of A. 
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1 
Let, v’ {;| . Using the iteration scheme (D.13), 


es) 


1 
sok, =2 and scaled v'= f 
0.5 


1 0 : a : 
lv —v |,,=0.5 > tol, so the iterations continue, 


fie 


Comparing the scaled vectors gives, 


lv’—v' |= 0.3 > tol so the iterations continue and after 6 iterations we have, 


1 
sok,=2.9918 and scaled v= 
0.9973 


7 7 6 3 7 . . 
Comparing scaled vectors gives, |v —v |,, = 0.0055 < tol, so the iterations stop. An estimate for the 


dominant eigenvalue of A is ks = 2.9918 (with v° being the corresponding estimated eigenvalue). The 
exact answer is 3. Computer code for this method can be downloaded from the website. 
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D.6 Rates of Convergence of Iterative Schemes 


The Rate of Convergence (RoC) of an iterative scheme is a measure of the number of iterations needed to 
converge to some given tolerance. It turns out that the RoC of an iterative scheme can be defined as, 


—log.% (D.14) 


where A is the dominant eigenvalue of its iteration matrix. 
The relative RoC of two schemes with dominant eigenvalues 1, and A, is, 


alog.A, (D.15) 
> log.A, 
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We can now compare (convergent) iterative schemes. 


e.g. Given that the point-Jacobi scheme is convergent with dominant eigenvalue A we know that the point- 
Gauss-Seidel scheme has dominant eigenvalue A’ and so is also convergent. By (D.15) the relative RoC of 
the Gauss-Seidel scheme to the Jacobi scheme is, 


2 
—log.A » 
—log.A 


i.e. the number of iterations to achieve the same level of accuracy using the Gauss-Seidel scheme is 
approximately half that of the Jacobi scheme. 
The relative RoC isn’t the whole story when comparing iterative schemes. It could be that an iterative 


scheme needs many iterations to converge but that each iteration is computationally fast. This could make 
it a faster than a quick converging but computationally slow scheme. 


D.7 Exercise D 


-—0.2 -03 -0.1 
1. 1. y_}_93 94. 935 |: Use Gershgorin’s theorem to provide an upper bound for the dominant 
-0.1 O05 0.5 


eigenvalue of A. If A was an iteration matrix would the iteration converge? 


21 


2. Find the exact dominant eigenvalue for 
12 


) and compare with the results from the power method 


code which you should run with a tolerance of 0.001. 


3. Adapt the power method code to estimate the dominant eigenvalue for the matrix in Q1. Check your 
answer by using Scilab’s (or Matlab’s) built-in function for finding eigenvalues. 


4. Given that the Jacobi iterative scheme converges show that the Gauss-Seidel scheme also converges. 
Show that the SoR scheme also converges. 


5. Find the RoC for an iterative scheme with dominant eigenvalue 0.3. 


6. Repeat Q5 for a dominant eigenvalue of 0.6. 
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7. Given that a Gauss-Seidel iterative scheme has dominant eigenvalue 0.5 find the relative RoC of SoR 
to Gauss-Seidel. 


8. Iterative scheme | has dominant eigenvalue 0.3 and iterative scheme 2 has dominant eigenvalue 0.6. 
Which scheme converges fastest? If scheme | takes 20 iterations to converge approximately in how 
many iterations does scheme 2 to converge (with the same tolerance)? 
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